Analisi di Clustering sul Dataset Banknote
D’Atri, Lombardi, Maiello
03/04/2025
1. Introduzione
Obiettivo del Progetto
Il presente progetto si propone di applicare diverse tecniche di clustering al dataset banknote al fine di identificare gruppi omogenei presenti nei dati e confrontare l’efficacia dei vari metodi. L’analisi intende mettere in pratica le metodologie apprese nel modulo di Data Analysis & Statistical Learning, in particolare, verranno trattati i seguenti algoritmi di clustering: Clustering Gerarchico Agglomerativo, Clustering Partizionale (K-means e K-medoids) e Clustering Model-based (Mixture of Gaussian).
TODO: specificare PCA e mixture of regression
Tra diversi dataset consigliati dai docenti, è stato scelto il dataset banknote, per un motivo puramente professionale, essendo tutti i membri del team legati professionalmente alle banconote. TODO: togliere?
Il documento è organizzato nei seguenti capitoli:
- Introduzione: Presentazione del progetto, degli obiettivi e del contesto.
- Descrizione del Dataset e Analisi Esplorativa: Esplorazione preliminare dei dati, analisi descrittiva e operazioni di preprocessing.
- Applicazione delle Tecniche di Clustering: Implementazione dei metodi di clustering (gerarchico, partizionale e model-based).
- Validazione e Visualizzazione dei Risultati: Confronto dei risultati mediante metriche di validazione e visualizzazione dei cluster (con l’ausilio della PCA).
- Conclusioni e Sviluppi Futuri: Sintesi dei risultati ottenuti e prospettive per ulteriori analisi.
2. Descrizione del Dataset e Analisi Esplorativa
In questa sezione esamineremo il dataset banknote
per comprenderne la struttura, le caratteristiche e le relazioni tra le
variabili. Il dataset proviene dal pacchetto mclust e
contiene misure che possono essere utilizzate per distinguere tra
banconote autentiche e false (o per altri scopi diagnostici). Il dataset
contiene le seguenti variabili:
Status: lo stato della banconota (autentica o contraffatta)
Length: lunghezza della banconota (mm)
Left: larghezza del bordo sinistro (mm)
Right: larghezza del bordo destro (mm)
Bottom: larghezza del margine inferiore (mm)
Top: larghezza del margine superiore (mm)
Diagonal: lunghezza della diagonale (mm)
Caricamento del Dataset
library(mclust)
data(banknote)Esplorazione della Struttura e Statistiche Descrittive del dataset
Visualizziamo le prime righe del dataset
head(banknote)## Status Length Left Right Bottom Top Diagonal
## 1 genuine 214.8 131.0 131.1 9.0 9.7 141.0
## 2 genuine 214.6 129.7 129.7 8.1 9.5 141.7
## 3 genuine 214.8 129.7 129.7 8.7 9.6 142.2
## 4 genuine 214.8 129.7 129.6 7.5 10.4 142.0
## 5 genuine 215.0 129.6 129.7 10.4 7.7 141.8
## 6 genuine 215.7 130.8 130.5 9.0 10.1 141.4Esaminiamo la struttura del dataset
str(banknote)## 'data.frame': 200 obs. of 7 variables:
## $ Status : Factor w/ 2 levels "counterfeit",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ Length : num 215 215 215 215 215 ...
## $ Left : num 131 130 130 130 130 ...
## $ Right : num 131 130 130 130 130 ...
## $ Bottom : num 9 8.1 8.7 7.5 10.4 9 7.9 7.2 8.2 9.2 ...
## $ Top : num 9.7 9.5 9.6 10.4 7.7 10.1 9.6 10.7 11 10 ...
## $ Diagonal: num 141 142 142 142 142 ...Il dataset è composto da 200 osservazioni e 7 variabili, tutte di tipo numerico ad eccezione di ‘Status’ che è una variabile binaria che può assumere uno tra i seguenti valori: ‘counterfeit’ e ‘genuine’.
Visualizziamo un riepilogo statistico delle variabili.
summary(banknote)## Status Length Left Right
## counterfeit:100 Min. :213.8 Min. :129.0 Min. :129.0
## genuine :100 1st Qu.:214.6 1st Qu.:129.9 1st Qu.:129.7
## Median :214.9 Median :130.2 Median :130.0
## Mean :214.9 Mean :130.1 Mean :130.0
## 3rd Qu.:215.1 3rd Qu.:130.4 3rd Qu.:130.2
## Max. :216.3 Max. :131.0 Max. :131.1
## Bottom Top Diagonal
## Min. : 7.200 Min. : 7.70 Min. :137.8
## 1st Qu.: 8.200 1st Qu.:10.10 1st Qu.:139.5
## Median : 9.100 Median :10.60 Median :140.4
## Mean : 9.418 Mean :10.65 Mean :140.5
## 3rd Qu.:10.600 3rd Qu.:11.20 3rd Qu.:141.5
## Max. :12.700 Max. :12.30 Max. :142.4Notiamo che la colonna ‘Status’ è perfettamente bilanciata.
Controlliamo se ci sono valori mancanti nelle colonne.
colSums(is.na(banknote))## Status Length Left Right Bottom Top Diagonal
## 0 0 0 0 0 0 0Come possiamo vedere nessuna colonna presenta valori mancanti.
Analisi Esplorativa dei Dati (EDA)
Questa sezione fornisce una panoramica iniziale del dataset banknote, evidenziando la struttura, le principali statistiche descrittive e le relazioni tra le variabili. Queste informazioni sono fondamentali per orientare le successive fasi di clustering e per capire eventuali necessità di preprocessing (ad es. scaling o trasformazioni) prima dell’applicazione dei metodi di clustering.
Analisi delle Correlazioni
Verifichiamo le correlazioni tra le variabili numeriche per capire se esistono relazioni forti che potrebbero influenzare l’analisi di clustering.
cor_matrix <- cor(banknote[, -1])
print(cor_matrix)## Length Left Right Bottom Top Diagonal
## Length 1.00000000 0.2312926 0.1517628 -0.1898009 -0.06132141 0.1943015
## Left 0.23129257 1.0000000 0.7432628 0.4137810 0.36234960 -0.5032290
## Right 0.15176280 0.7432628 1.0000000 0.4867577 0.40067021 -0.5164755
## Bottom -0.18980092 0.4137810 0.4867577 1.0000000 0.14185134 -0.6229827
## Top -0.06132141 0.3623496 0.4006702 0.1418513 1.00000000 -0.5940446
## Diagonal 0.19430146 -0.5032290 -0.5164755 -0.6229827 -0.59404464 1.0000000library(reshape2)
library(ggplot2)
cor_melted <- melt(cor_matrix)
ggplot(cor_melted, aes(x = Var1, y = Var2, fill = value)) + geom_tile() +
geom_text(aes(label = round(value, 2)), color = "black",
size = 5) + scale_fill_gradient2(low = "blue", mid = "white",
high = "red", midpoint = 0, limits = c(-1, 1)) + theme_minimal() +
labs(title = "Correlazione variabili", fill = "Correlazione") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
TODO: Commento:
Osserva i valori di correlazione per identificare variabili fortemente correlate (positivamente o negativamente).
Discuti l’impatto di eventuali correlazioni elevate: potrebbero ridondare in informazioni e influenzare il clustering.
Se trovi correlazioni elevate, suggerisci che una riduzione della dimensionalità (es. tramite PCA) potrebbe aiutare.Distribuzione delle Variabili
Infine, esaminiamo la distribuzione di ciascuna variabile con degli istogrammi:
# Impostiamo un layout per visualizzare gli istogrammi
par(mfrow = c(2, 2))
for(i in 2:7){
hist(banknote[[i]], main = paste("Istogramma di", names(banknote)[i]),
xlab = names(banknote)[i], col = "lightblue", border = "white")
}
par(mfrow = c(1, 1)) # Ripristiniamo il layout originale
TODO: Commento:
Descrivi la forma delle distribuzioni: simmetriche, asimmetriche (skewed), unimodali o multimodali.
Nota eventuali outlier o code lunghe che potrebbero influenzare la successiva analisi di clustering.
Confronta le distribuzioni tra le variabili per capire se alcune variabili hanno una variabilità molto diversa, che potrebbe richiedere una normalizzazione/scaling.Matrice di Scatterplot
Per avere una prima idea delle relazioni tra le variabili, creiamo una matrice di scatterplot. Utilizziamo la colonna ‘Status’ per colorare i punti.
library(GGally)
ggpairs(banknote[,-1],upper = list(continuous = "density", combo = "box_no_facet"),
lower = list(continuous = "points", combo = "dot_no_facet"),aes(colour=banknote$Status))
Lungo la diagonale sono visibili le distribuzioni univariate delle variabili. Si nota una sovrapposizione parziale tra i due gruppi di colore (Status), ma anche zone in cui uno dei due colori prevale. Nei grafici di densità bivariata (metà diagonale superiore della matrice) e negli scatter plot (metà diagonale inferiore della matrice), si osservano alcune regioni in cui i punti di uno stesso colore tendono a concentrarsi, suggerendo cluster parzialmente distinti. Nonostante la sovrapposizione in diverse aree, si intravedono regioni in cui un gruppo prevale. Ciò suggerisce che un algoritmo di classificazione o di clustering potrebbe distinguere parzialmente le due classi, anche se non in modo perfetto con una singola coppia di variabili. È importante sottolineare che questo grafico mostra solo le relazioni a coppie. Ciò significa che l’analisi si limita a considerare due variabili alla volta, mentre una clusterizzazione netta potrebbe emergere solo quando si considerano contemporaneamente più variabili. In altre parole, anche se alcune coppie non evidenziano una separazione chiara, combinando più dimensioni si potrebbe scoprire una struttura clusterizzata più marcata. Il fatto che in alcune proiezioni (coppie di variabili) si vedano zone di separazione o contorni distinti suggerisce che esistono sottostrutture nei dati. Notiamo infine che le variabili ‘Bottom’ e ‘Diagonal’ sembrerebbero separare in modo migliore il dataset rispetto alle altre variabili.
D’ora in avanti, per ciascun algoritmo di clustering che proveremo, verranno effettuate tre diverse analisi:
Utilizzando tutte le variabili disponibili nel dataset.
Considerando esclusivamente le variabili ‘Diagonal’ e ‘Bottom’.
Applicando l’Analisi delle Componenti Principali (PCA) per ridurre la dimensionalità.
Queste tre approcci permetteranno di valutare l’impatto delle diverse selezioni di variabili e della riduzione dimensionale sulle prestazioni e sui risultati degli algoritmi di clustering.
Prima di applicare gli algoritmi di clustering, è fondamentale verificare se il dataset presenta una naturale propensione alla formazione di gruppi.
Valutazione della tendenza al clustering del dataset banknote
Utilizziamo specifiche misure statistiche per determinare se i dati sono strutturati in cluster o se, al contrario, sono distribuiti in modo casuale. Questi metodi permettono di valutare se il clustering è effettivamente significativo e utile per l’analisi del dataset.
Hopkins statistics
L’Hopkins Statistic è una misura utilizzata per valutare la tendenza di un dataset alla clusterizzazione, confrontando la sua struttura con quella di un insieme di punti generati casualmente. Assume valori compresi tra 0 e 1. Valori vicini a 0 indicano che il dataset è altamente strutturato in cluster, suggerendo che l’applicazione di algoritmi di clustering può essere efficace. Valori vicini a 0.5 o superiori suggeriscono che i dati sono distribuiti in modo casuale o uniforme, e il clustering potrebbe non essere significativo.
library(clustertend)
banknote_num <- banknote[, -1]
banknote_num <- scale(banknote_num)
set.seed(987)
hopkins_stat <- hopkins(banknote_num, n = nrow(banknote_num) - 1)
random_df <- apply(banknote_num, 2, function(x){runif(length(x), min(x), max(x))})
random_df <- as.data.frame(random_df)
random_df <- scale(random_df)
set.seed(987)
hopkins_stat_rand <- hopkins(random_df, n = nrow(random_df)-1)
cat("Hopkins Statistics on Banknote dataset:", hopkins_stat$H, "\n")## Hopkins Statistics on Banknote dataset: 0.2642237cat("Hopkins Statistics on Random dataset:", hopkins_stat_rand$H, "\n")## Hopkins Statistics on Random dataset: 0.4926061Un valore di Hopkins pari a 0.2519264, essendo inferiore a 0.5, indica una tendenza alla clusterizzazione. In altre parole, le distanze tra i punti reali sono molto più ridotte rispetto a quelle dei punti generati casualmente, suggerendo la presenza di gruppi ben definiti nel dataset. Questo risultato supporta l’ipotesi che il dataset sia strutturato in cluster distinti.
Visual Assessment of cluster Tendency (VAT)
L’algoritmo VAT produce una matrice di dissimilarità ordinata, ottenuta riorganizzando le righe e le colonne della matrice delle distanze per evidenziare eventuali blocchi di punti “ravvicinati” (ossia potenziali cluster). Visualizziamo la matrice con l’intento di verificare se si formano lungo la diagonale secondaria del grafico, dei blocchi rossi che rappresentano regioni di osservazioni che condividono una forte similarità.
library("factoextra")
library(gridExtra)
p1 <- fviz_dist(dist(banknote), show_labels = FALSE) + labs(title = "Banknote data")
p2 <- fviz_dist(dist(random_df), show_labels = FALSE) + labs(title = "Random data")
grid.arrange(p1, p2, ncol = 2, widths = c(6, 6))
Nei grafici VAT il rosso indica alta similarità (ossia bassa
dissimilarità), mentre il blu indica bassa similarità (ossia alta
dissimilarità). Relativamente al grafico del dataset banknote, notiamo
che lungo la diagonale sono presenti regioni rosse che indicano dati
simili tra loro, suggerendo la presenza di cluster. La presenza di zone
blu che evidenziano una maggiore dissimilarità tra gruppi di dati,
supporta l’ipotesi che esistano cluster distinti nel dataset. Inoltre la
presenza di transizioni marcate dal rosso al blu conferma che esistono
confini ben definiti tra i cluster. Notiamo invece come nel dataset
randomico non c’è alcuna struttura a cluster evidente. I dati appaiono
disordinati e non presenteranno raggruppamenti chiari, il che indica che
un clustering non ha senso o che i dati non sono adatti per essere
separati in gruppi distinti.
Possiamo concludere che il dataset Banknote ha una tendezza alla clusterizzazione, procediamo perciò con il calcolo della PCA e successivamente alla clusterizzazione.
PCA
In questa sezione applichiamo la PCA sul dataset, considerando solo le colonne numeriche. Anche se le variabili sono tutte misurate in millimetri e quindi nella stessa scala, applichiamo comunque la standardizzazione per evitare che differenze nella dispersione (variabili con varianze maggiori) influenzino in modo sproporzionato i risultati della PCA.
library("FactoMineR")
banknote_num <- banknote[, -1]
head(banknote_num)## Length Left Right Bottom Top Diagonal
## 1 214.8 131.0 131.1 9.0 9.7 141.0
## 2 214.6 129.7 129.7 8.1 9.5 141.7
## 3 214.8 129.7 129.7 8.7 9.6 142.2
## 4 214.8 129.7 129.6 7.5 10.4 142.0
## 5 215.0 129.6 129.7 10.4 7.7 141.8
## 6 215.7 130.8 130.5 9.0 10.1 141.4pca_result <- PCA(banknote_num, scale = TRUE, graph = FALSE)Visualizzazione e interpretazione dei risultati
Estraiamo gli autovalori e calcoliamo sia la varianza spiegata per ciascuna componente sia la varianza cumulativa. Queste informazioni ci serviranno per decidere il numero ottimale di componenti principali da mantenere.
eigenvalues <- get_eigenvalue(pca_result)[,1]
explained_variance <- eigenvalues / sum(eigenvalues)
cumulative_variance <- cumsum(explained_variance)
print(eigenvalues)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5 Dim.6
## 2.9455582 1.2780838 0.8690326 0.4497687 0.2686769 0.1888799print(explained_variance)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5 Dim.6
## 0.49092637 0.21301396 0.14483876 0.07496145 0.04477948 0.03147998print(cumulative_variance)## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5 Dim.6
## 0.4909264 0.7039403 0.8487791 0.9237405 0.9685200 1.0000000Abbiamo ora a disposizione gli autovalori e le informazioni sulla varianza, utili per valutare l’importanza di ciascuna componente principale.
Utilizziamo la Kaiser Rule, una tecnica utilizzata per determinare il numero di componenti principali da mantenere quando si esegue l’analisi delle componenti principali. La regola si basa sull’idea di conservare solo le componenti con una varianza spiegata maggiore di 1. Questo criterio si applica al valore proprio (eigenvalue) delle componenti principali.
kaiser_components <- sum(eigenvalues > 1)
cat("Numero di componenti selezionati con la Kaiser rule:", kaiser_components, "\n")## Numero di componenti selezionati con la Kaiser rule: 2Otteniamo 2 come numero di componenti principali da mantenere, ciò significa che le due componenti principali con valori propri superiori a 1 spiegano una quantità di varianza maggiore di quella che verrebbe spiegata da una singola variabile originale del dataset.
Un ulteriore criterio è quello di selezionare il numero minimo di componenti la cui somma di varianza spiegata (varianza comulativa), raggiunge almeno l’80% del totale. In questo modo garantiamo una buona rappresentazione dei dati.
variance_threshold <- min(which(cumulative_variance >= 0.80))
cat("Numero di componenti principali che spiegano l'80% della varianza comulativa:", variance_threshold, "\n")## Numero di componenti principali che spiegano l'80% della varianza comulativa: 3fviz_eig(pca_result, addlabels = TRUE, ylim = c(0, 50))
Nel grafico vediamo la varianza spiegata da ciascuna componente. Con le
prime tre componenti principali riusciamo a spiegare l’84.9% della
varianza cumulativa.
Vediamo infine un ultimo metodo per scegliere il numero di componenti principali: l’elbow method. Mediante un grafico scree plot visualizziamo gli autovalori per ciascuna componente e individuiamo “il gomito”, ovvero il punto nel grafico in cui la diminuzione degli autovalori rallenta. Tale punto può essere interpretato come un’indicazione del numero ottimale di componenti. Includiamo nel grafico anche due linee verticali per evidenziare i criteri della varianza cumulativa (rosso) e della Kaiser Rule (blu).
scree_plot <- data.frame(PC = 1:length(eigenvalues), Eigenvalue = eigenvalues)
ggplot(scree_plot, aes(x = PC, y = Eigenvalue)) +
geom_point() +
geom_line() +
geom_vline(xintercept = variance_threshold, linetype = "dashed", color = "red") +
geom_vline(xintercept = kaiser_components, linetype = "dotted", color = "blue") +
labs(title = "Scree Plot for PCA", x = "Principal Component", y = "Eigenvalue") +
theme_minimal()
Nel grafico il gomito sembra collocarsi intorno alla quarta
componente. I tre criteri danno perciò risultati differenti, per
mantenere un equilibrio tra interpretabilità e preservazione della
varianza, si potrebbe optare per 3 componenti.
Di seguito viene costruito il nuovo dataset in cui le variabili sono
sostituite con le prime 3 PCA, che verrà usato più avanti nelle
analisi.
pca_data <- as.data.frame(pca_result$ind$coord[, 1:3])
colnames(pca_data) <- c("PCA1", "PCA2", "PCA3")
banknote_pca <- cbind(pca_data, Status = banknote$Status)
banknote_pca <- banknote_pca %>% dplyr::select(Status, everything())
head(banknote_pca)## Status PCA1 PCA2 PCA3
## 1 genuine 1.7474012 1.65082829 1.4237611
## 2 genuine -2.2743177 -0.53879328 0.5326484
## 3 genuine -2.2774015 -0.10767707 0.7174149
## 4 genuine -2.2835546 -0.08765431 -0.6056336
## 5 genuine -2.6321283 0.03919590 3.1963847
## 6 genuine 0.7584073 3.08874513 0.7864804Per comprendere meglio il contributo delle variabili alle componenti principali, visualizziamo le loro coordinate e realizziamo una mappa di correlazione. Questo ci aiuta a capire come le variabili si relazionano alle componenti principali.
var <- get_pca_var(pca_result)
var$coord## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## Length -0.01199158 0.9219364 -0.01648225 -0.38536590 0.0304764
## Left 0.80279596 0.3866019 0.09637548 0.26485400 -0.3314768
## Right 0.83526859 0.2854104 0.11510550 0.28856524 0.3183114
## Bottom 0.69810421 -0.3009779 0.54398559 -0.27072283 0.1116898
## Top 0.63139786 -0.1034278 -0.73418921 -0.07392333 0.1139569
## Diagonal -0.84690418 0.3096965 0.10615680 0.26284740 0.1763188fviz_pca_var(pca_result, col.var = "contrib",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"))
Questo passaggio ci offre una visione immediata di come ogni variabile
contribuisce alle componenti principali e della loro interrelazione.
Approfondiamo l’analisi mostrando i contributi delle variabili attraverso diversi metodi grafici, in modo da identificare quali variabili influenzano maggiormente ciascuna componente.
print(round(var$contrib,digits = 3))## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## Length 0.005 66.503 0.031 33.019 0.346
## Left 21.880 11.694 1.069 15.596 40.896
## Right 23.686 6.374 1.525 18.514 37.712
## Bottom 16.545 7.088 34.052 16.295 4.643
## Top 13.534 0.837 62.027 1.215 4.833
## Diagonal 24.350 7.504 1.297 15.361 11.571library("corrplot")
corrplot(var$contrib, is.corr=FALSE)
L’analisi dei contributi evidenzia il peso relativo di ciascuna
variabile nelle componenti principali, utile per interpretare i
risultati della PCA.
Visualizziamo i contributi delle variabili specificamente per la prima componente principale.
fviz_contrib(pca_result, choice = "var", axes = 1)
Similmente al passaggio precedente, visualizziamo i contributi delle
variabili per la seconda componente principale.
fviz_contrib(pca_result, choice = "var", axes = 2)
3. Applicazione delle Tecniche di Clustering
3.1 Clustering Gerarchico Agglomerativo
Il Clustering Gerarchico Agglomerativo (HAC) è una tecnica di clustering che costruisce una gerarchia di cluster attraverso un approccio “bottom-up”. Inizia considerando ogni punto dati come un cluster separato e, successivamente, unisce iterativamente i cluster più simili fino a formare un unico cluster che racchiude tutti i dati. Questo metodo non richiede la specifica del numero di cluster a priori e produce una rappresentazione ad albero, nota come dendrogramma, che facilita la visualizzazione delle relazioni tra i cluster.
Per applicare il clustering gerarchico agglomerativo al nostro dataset, seguiamo questi passaggi:
Scaling:
Calcolo della Matrice di Distanza: Determiniamo le distanze tra tutte le coppie di punti nel dataset. Una misura comune è la distanza euclidea.
Applicazione del clustering agglomerativo: Utilizziamo la funzione
hclust()per eseguire il clustering gerarchico e generare il dendrogramma.Taglio del Dendrogramma: Decidiamo il numero ottimale di cluster utilizzando varie tecniche, e tagliamo il dendrogramma al livello appropriato.
Valutazione: Valutiamo i risultati del clustering utilizzando varie metriche.
Clustering con tutte le variabili
Scaliamo il dataset e calcoliamo la relativa matrice di dissimilatità mediante distanza euclidea e la mostriamo. si è scelta la distanza euclidea poichè tutte le variabili del dataset su cui applicheremo il clustering sono di tipo numerico.
library(pheatmap)
df <- scale(banknote_num)
distanze <- dist(df, method = "euclidean")
diss_matrix <- as.matrix(distanze)
pheatmap(diss_matrix,
cluster_rows = FALSE,
cluster_cols = FALSE,
show_rownames = FALSE,
show_colnames = FALSE,
color = colorRampPalette(c("blue", "white", "red"))(50),
border_color = NA)
Applichiamo il clustering con metodo di linkage ‘ward’ che generalmente insieme all’average linkage, performa meglio.
hclust_clustering <- hclust(distanze, method = "ward.D2")Guardiamo il dendogramma cercando di identificare il numero ottimale di cluster. Utilizzeremo poi anche altri metodi che ci premetteranno di prendere una scelta migliore.
Numero di cluster
fviz_dend(hclust_clustering, show_labels = FALSE, as.ggplot = TRUE)
Il dendrogramma risultante mostra come le osservazioni vengono unite
progressivamente in cluster. L’altezza dei rami indica la distanza a cui
avviene la fusione: rami più bassi corrispondono a fusioni tra punti o
cluster molto simili. Analizzando il dendrogramma, possiamo decidere il
numero di cluster tagliando l’albero a un’altezza che separa i gruppi in
modo significativo. Dall’osservazione del dendrogramma, si nota un salto
particolarmente ampio a un’altezza intorno a 30, il che indica che fino
a quella soglia i gruppi sono relativamente compatti e ben separati.
Tagliare il dendrogramma a quel livello produrrebbe 2 grandi cluster.
Volendo una suddivisione più fine, si può individuare un altro salto
meno marcato (intorno a 15) che porterebbe a 3 cluster. Procediamo con
l’analisi del numero ottimale di cluster utilizzando altri metodi.
Osserviamo l’andamento del Total Within Sum of Squares (WSS) rispetto a k, cercando di individuare il punto in cui la curva inizia a “piegarsi” più marcatamente, cioè dove l’aggiunta di ulteriori cluster non porta a una riduzione significativa del WSS. Di seguito il grafico Elbow.
fviz_nbclust(df, FUN = hcut, method = "wss") +
geom_vline(xintercept = 3, linetype = 2) +
labs(subtitle = "Elbow method")
Tra k = 1 e k = 2, la diminuzione del WSS è molto forte (ciò è normale,
perché si passa da un unico cluster a due). Tra k = 2 e k = 3, c’è
ancora un calo netto del WSS. Oltre k = 3, la pendenza inizia a ridursi
gradualmente. Dal grafico, sembra che il calo principale avvenga per k =
3, mentre oltre si riduce più lentamente. Per k = 3 si ha un notevole
miglioramento rispetto a 2, ma passando a 4 (e oltre) il guadagno in
termini di riduzione del WSS diventa più contenuto.
Analizziamo ora il grafico dell’Average Silhouette Width.
fviz_nbclust(df, FUN = hcut, method = "silhouette") +
labs(subtitle = "Silhouette method")
L’indice di silhouette risulta più alto in corrispondenza di k = 2.
Questo suggerisce che la separazione tra i cluster e la coesione interna
sia massima quando si dividono i dati in due gruppi. Con k = 3, l’indice
diminuisce, e oltre tale valore rimane su livelli inferiori (intorno a
0.2-0.3), senza tornare ai valori di picco.
Osserviamo il grafico del Gap Statistic che confronta la dispersione all’interno dei cluster ottenuta dal clustering, con quella attesa in un dataset generato casualmente (quindi privo di gruppi). Ossia calcola la differenza (gap) tra il log della dispersione osservata e quella media attesa sotto una distribuzione nulla.
set.seed(123)
fviz_nbclust(df, FUN = hcut, method = "gap_stat", nboot = 500)+
labs(subtitle = "Gap statistic method")
Nel grafico, la linea che rappresenta il punto in cui il Gap Statistic
non aumenta più in maniera significativa, è su 3 cluster, il che indica
che quando i dati vengono divisi in 3 gruppi la struttura di clustering
risulta migliore rispetto a quella attesa in un dataset casuale. In
altre parole, il modello a 3 cluster spiega meglio la struttura
intrinseca dei dati.
Utilizziamo infine il pacchetto ‘NbClust’ che valuta il numero ottimale di cluster utilizzando numerosi indici interni. Tra questi, alcuni indici come l’indice di Hubert e il Dindex vengono visualizzati graficamente.
library("NbClust")
nb <- NbClust(df, distance = "euclidean", min.nc = 2, max.nc= 10, method = "ward.D2")
## *** : The Hubert index is a graphical method of determining the number of clusters.
## In the plot of Hubert index, we seek a significant knee that corresponds to a
## significant increase of the value of the measure i.e the significant peak in Hubert
## index second differences plot.
## 
## *** : The D index is a graphical method of determining the number of clusters.
## In the plot of D index, we seek a significant knee (the significant peak in Dindex
## second differences plot) that corresponds to a significant increase of the value of
## the measure.
##
## *******************************************************************
## * Among all indices:
## * 7 proposed 2 as the best number of clusters
## * 10 proposed 3 as the best number of clusters
## * 1 proposed 4 as the best number of clusters
## * 1 proposed 7 as the best number of clusters
## * 1 proposed 8 as the best number of clusters
## * 3 proposed 10 as the best number of clusters
##
## ***** Conclusion *****
##
## * According to the majority rule, the best number of clusters is 3
##
##
## *******************************************************************La decisione finale basata sulla majority rule suggerisce 3 cluster, indicando che, complessivamente, questa soluzione offre una migliore separazione e coesione interna rispetto alle altre soluzioni.
Complessivamente la maggioranza dei metodi analizzati (elbow, gap statistic, nbclust) converge su una soluzione a 3 cluster. Pertanto optiamo per una scelta di “maggioranza” e scegliamo perciò un numero di cluster pari a 3.
Taglio del dendrogramma
Procediamo con il taglio del dendogramma in modo da ottenere 3 cluster. Successivamente visualizziamo una tabella di frequenza che mostra il numero di osservazioni assegnate a ciascun cluster, permettendoci di verificare la distribuzione dei dati tra i 3 cluster ottenuti.
cluster_cut <- cutree(hclust_clustering, k = 3)
table(cluster_cut)## cluster_cut
## 1 2 3
## 102 76 22fviz_dend(hclust_clustering, k = 3,
cex = 0.5,
k_colors = c("#2E9FDF", "#00AFBB", "#E7B800", "#FC4E07"),
color_labels_by_k = TRUE,
rect = TRUE
)
La maggior parte dei dati è stata assegnata ai primi due cluster, mentre
il terzo cluster contiene un numero significativamente minore di
osservazioni. Questo potrebbe suggerire che il cluster 3 rappresenta un
gruppo più piccolo, con caratteristiche potenzialmente diverse o più
specifiche rispetto ai gruppi più grandi.
Vediamo i risultati di clustering nello spazio originale, tramite la matrice di scatterplot a coppie con i punti colorati in base ai cluster di appartenenza.
ggpairs(df, aes(color = as.factor(cluster_cut))) +
ggtitle("Matrice Pairwise Scatterplot con Cluster")
Valutazione
Valutiamo ora la qualità dei risultati del nostro algoritmo di clustering mediante quattro metriche: Average Silhouette Width, Dunn Index, Corrected Rand Index e Meila’s Variation of Information Index.
Iniziamo con la metrica della silhouette media (calcolata come la media di s(i) su tutte le osservazioni), che fornisce un’indicazione complessiva della qualità del clustering: si auspica un valore il più vicino ad 1 possibile.
library(cluster)
sil <- silhouette(cluster_cut, distanze)
fviz_silhouette(sil, palette = "jco", ggtheme = theme_classic())## cluster size ave.sil.width
## 1 1 102 0.30
## 2 2 76 0.30
## 3 3 22 0.39
Un valore di 0.31 suggerisce che la struttura dei cluster non è
particolarmente forte: c’è una separazione moderata, ma c’è anche una
certa sovrapposizione tra i cluster.
Analizziamo ora il Dunn Index, che valuta la separazione tra i cluster rispetto alla loro compattezza, cioè quanto i cluster siano ben separati e internamente compatti. Se i dati contengono cluster compatti e ben separati, ci si aspetta che il diametro dei cluster sia piccolo e che la distanza tra i cluster sia grande. Pertanto l’indice di Dunn deve essere massimizzato.
library(clusterCrit)
dunn_index <- intCriteria(as.matrix(df), cluster_cut, "Dunn")$dunn
print(paste("Dunn Index:", dunn_index))## [1] "Dunn Index: 0.1250576993898"Un Dunn Index di 0.125 indica che i cluster hanno una separazione piuttosto bassa rispetto alla loro coesione interna. Una bassa separazione tra cluster suggerisce possibili problemi di sovrapposizione.
Analizziamo il Corrected Rand Index che misura quanto bene il clustering ottenuto corrisponde a una partizione di riferimento. Nel nostro caso utilizziamo la colonna ‘Status’ del dataset. Ed infine il Meila’s Variation of Information (VI) Index che misura quanto due cluster differiscono tra loro in termini di entropia. A differenza del Corrected Rand Index, che valuta la similarità, il VI Index misura la distanza tra due partizioni (va perciò minimizzato).
library("fpc")
status <- as.numeric(banknote$Status)
clust_stats <- cluster.stats(d = dist(df), status, cluster_cut)
print(paste("Corrected Rand Index:", clust_stats$corrected.rand))## [1] "Corrected Rand Index: 0.791980314729782"print(paste("Meila's Index:", clust_stats$vi))## [1] "Meila's Index: 0.359179851985559"Un Corrected Rand Index di 0.792 indica che il clustering è molto simile alla vera classificazione dei dati. Un Meila’s Variation of Information Index di 0.359 è un valore relativamente basso, il che significa che i cluster trovati non sono troppo distanti dalla struttura effettiva dei dati. Il fatto che il Corrected Rand Index sia alto (0.792) e il VI Index sia basso conferma che il clustering ha una buona coerenza con la partizione reale. Tuttavia, poiché Dunn Index e Silhouette suggeriscono che la separazione tra cluster è debole, questo potrebbe indicare sovrapposizione tra i cluster, ma senza compromettere eccessivamente la qualità dell’assegnazione.
Ottimizzazione
Ottimizziamo il clustering esplorando diverse combinazioni di metodi di linkage e numeri di cluster. Per ciascuna combinazione, calcoliamo le quattro metriche di valutazione precedentemente analizzate, ottenendo un dataframe con tutti i risultati. Successivamente, identifichiamo il metodo di clustering migliore per ciascuna metrica. Poiché le diverse metriche potrebbero suggerire metodi differenti, assegniamo un ranking a ciascuna di esse e calcoliamo un ranking totale sommando i punteggi (invertendo il criterio per le metriche da massimizzare). Infine, ordiniamo i risultati e selezioniamo la combinazione con il punteggio complessivo più basso, individuando così il miglior clustering in base a tutte le metriche disponibili.
library(dplyr)
evaluate_clustering <- function(df, distanze, status, linkage_methods = c("ward.D2", "single", "complete", "average", "centroid"), min_k = 2, max_k = 10) {
results_all <- data.frame(Method = character(),
k = numeric(),
Silhouette = numeric(),
Dunn = numeric(),
Corrected_Rand = numeric(),
VI = numeric(),
stringsAsFactors = FALSE)
for(method in linkage_methods) {
hclust_obj <- hclust(distanze, method = method)
for(k in min_k:max_k) {
cluster_cut <- cutree(hclust_obj, k = k)
sil <- silhouette(cluster_cut, distanze)
avg_sil <- mean(sil[, 3])
dunn_index <- intCriteria(as.matrix(df), cluster_cut, "Dunn")$dunn
clust_stats <- cluster.stats(distanze, status, cluster_cut)
corrected_rand <- clust_stats$corrected.rand
vi_val <- clust_stats$vi
results_all <- rbind(results_all, data.frame(Method = method,
k = k,
Silhouette = avg_sil,
Dunn = dunn_index,
Corrected_Rand = corrected_rand,
VI = vi_val,
stringsAsFactors = FALSE))
}
}
return(results_all)
}
results_all <- evaluate_clustering(df, distanze, status)
best_silhouette <- results_all %>% group_by(Method) %>% filter(Silhouette == max(Silhouette))
print(best_silhouette)## # A tibble: 5 × 6
## # Groups: Method [5]
## Method k Silhouette Dunn Corrected_Rand VI
## <chr> <int> <dbl> <dbl> <dbl> <dbl>
## 1 ward.D2 2 0.375 0.211 0.960 0.0982
## 2 single 2 0.333 0.341 0 0.718
## 3 complete 2 0.346 0.125 0.591 0.566
## 4 average 2 0.333 0.341 0 0.718
## 5 centroid 2 0.333 0.341 0 0.718best_dunn <- results_all %>% group_by(Method) %>% filter(Dunn == max(Dunn))
print(best_dunn)## # A tibble: 5 × 6
## # Groups: Method [5]
## Method k Silhouette Dunn Corrected_Rand VI
## <chr> <int> <dbl> <dbl> <dbl> <dbl>
## 1 ward.D2 2 0.375 0.211 0.960 0.0982
## 2 single 2 0.333 0.341 0 0.718
## 3 complete 10 0.186 0.179 0.268 1.45
## 4 average 2 0.333 0.341 0 0.718
## 5 centroid 2 0.333 0.341 0 0.718best_corrected_rand <- results_all %>% group_by(Method) %>% filter(Corrected_Rand == max(Corrected_Rand))
print(best_corrected_rand)## # A tibble: 5 × 6
## # Groups: Method [5]
## Method k Silhouette Dunn Corrected_Rand VI
## <chr> <int> <dbl> <dbl> <dbl> <dbl>
## 1 ward.D2 2 0.375 0.211 0.960 0.0982
## 2 single 8 -0.0124 0.225 0.000320 0.896
## 3 complete 3 0.327 0.138 0.688 0.596
## 4 average 5 0.316 0.222 0.951 0.140
## 5 centroid 3 0.228 0.298 0.000101 0.742best_vi <- results_all %>% group_by(Method) %>% filter(VI == min(VI))
print(best_vi)## # A tibble: 5 × 6
## # Groups: Method [5]
## Method k Silhouette Dunn Corrected_Rand VI
## <chr> <int> <dbl> <dbl> <dbl> <dbl>
## 1 ward.D2 2 0.375 0.211 0.960 0.0982
## 2 single 2 0.333 0.341 0 0.718
## 3 complete 2 0.346 0.125 0.591 0.566
## 4 average 5 0.316 0.222 0.951 0.140
## 5 centroid 2 0.333 0.341 0 0.718library(dplyr)
results_ranked <- results_all %>%
mutate(
rank_sil = rank(-Silhouette, ties.method = "min"),
rank_dunn = rank(-Dunn, ties.method = "min"),
rank_corrected_rand = rank(-Corrected_Rand, ties.method = "min"),
rank_vi = rank(VI, ties.method = "min"),
total_rank = rank_sil + rank_dunn + rank_corrected_rand + rank_vi
)
best_overall <- results_ranked %>% arrange(total_rank) %>% head(1)
print(best_overall)## Method k Silhouette Dunn Corrected_Rand VI rank_sil rank_dunn
## 1 ward.D2 2 0.3748428 0.2107915 0.9602001 0.09823913 1 21
## rank_corrected_rand rank_vi total_rank
## 1 1 1 24Il risultato mostra che il clustering ottimale secondo queste metriche è Ward.D2 con 2 cluster, poiché bilancia al meglio la coesione, la separazione e la corrispondenza con la verità di riferimento. La combinazione scelta restituisce i valori migliori possibili per le metriche Silhouette, Corrected Rand Index e VI (Rank 1), mentre risulta meno performante rispetto ad altri metodi, ma compensato dagli altri punteggi, per il Dunn Index (Rank 21).
Effettuiamo il clustering con la combinazione risultante e visualizziamo mediante il grafico i cluster ottenuti.
final <- eclust(df, "hclust", k = 2, hc_metric = "euclidean", hc_method = "ward.D2", graph = FALSE)
fviz_cluster(final, geom = "point", ellipse.type = "norm", palette="jco", ggtheme = theme_minimal())
Clustering con PCA
Vediamo come cambia il clustering se invece che considerare tutte le variabili del dataset, consideriamo solo le tre componenti principali trovate prima.
banknote_pca_num <- banknote_pca[,-1]
results_all_pca <- evaluate_clustering(banknote_pca_num, dist(banknote_pca_num), status)
results_ranked_pca <- results_all_pca %>%
mutate(
rank_sil = rank(-Silhouette, ties.method = "min"),
rank_dunn = rank(-Dunn, ties.method = "min"),
rank_corrected_rand = rank(-Corrected_Rand, ties.method = "min"),
rank_vi = rank(VI, ties.method = "min"),
total_rank = rank_sil + rank_dunn + rank_corrected_rand + rank_vi
)
best_overall_pca <- results_ranked_pca %>% arrange(total_rank) %>% head(1)
print(best_overall_pca)## Method k Silhouette Dunn Corrected_Rand VI rank_sil rank_dunn
## 1 average 2 0.431344 0.1056305 0.9020088 0.1997658 1 27
## rank_corrected_rand rank_vi total_rank
## 1 2 1 31Il risultato indica che, applicando il clustering agglomerativo con numero di cluster pari a 2 e metodo di linkage “average” al dataset ridotto a tre componenti principali, otteniamo la combinazione migliore in base ai 4 indici di valutazione, anche se i risultati presentano alcune peculiarità. Rispetto alle altre combinazioni testate, la silhouette e la VI hanno ottenuto i migliori punteggi (rank 1), il Corrected Rand è quasi il migliore (rank 2), invece l’indice Dunn, in questa combinazione, non performa bene rispetto ad altre configurazioni testate (rank 27). Rispetto al clustering agglomerativo applicato su tutte le variabili del dataset il clustering sul dataset ridotto mostra: - Un indice silhouette leggermente migliore (0.431) rispetto a quello su tutte le variabili (0.375), suggerendo che i punti sono, mediamente, meglio raggruppati all’interno dei cluster nel caso PCA. - Il clustering su tutte le variabili ha un indice Dunn migliore (0.211 vs 0.106), il che significa che la separazione minima tra cluster è maggiore e i cluster risultano più distinti. - Entrambe le soluzioni mostrano valori molto alti del Corrected Rand e bassi del VI, ma la soluzione con tutte le variabili è leggermente superiore (0.960 vs 0.902 e 0.098 vs 0.200), indicando una migliore aderenza ad una struttura di riferimento o una maggiore stabilità del clustering.
Mentre l’approccio con PCA (riduzione a 3 componenti) offre cluster con una buona compattezza (indice silhouette migliore), il clustering basato su tutte le variabili, utilizzando il metodo ward.D2, risulta complessivamente superiore per quanto riguarda la separazione tra cluster e la corrispondenza con una struttura di riferimento (migliori valori di Dunn, Corrected Rand e VI).
Vediamo graficamente come sono separati i due cluster.
library(scatterplot3d)
final_pca <- eclust(banknote_pca_num, "hclust", k = 2, hc_metric = "euclidean", hc_method = "average", graph = FALSE)
clusters <- final_pca$cluster
scatterplot3d(banknote_pca_num,
color = clusters,
pch = 19,
main = "Clustering su PCA (PC1, PC2, PC3)",
xlab = "PC1", ylab = "PC2", zlab = "PC3")
legend("topright", legend = unique(clusters), col = unique(clusters), pch = 19, title = "Cluster")
Clustering con solo le variabili Bottom e Diagonal
Vediamo infine come cambia il clustering considerando solo le variabili ‘Bottom’ e ‘Diagonal’ che sembravano suddividere meglio i dati in gruppi.
banknote_bottom_diagonal <- banknote %>% select(Bottom, Diagonal)
banknote_bd_scaled <- scale(banknote_bottom_diagonal)
results_all_bd <- evaluate_clustering(banknote_bd_scaled, dist(banknote_bd_scaled), status)
results_ranked_bd <- results_all_bd %>%
mutate(
rank_sil = rank(-Silhouette, ties.method = "min"),
rank_dunn = rank(-Dunn, ties.method = "min"),
rank_corrected_rand = rank(-Corrected_Rand, ties.method = "min"),
rank_vi = rank(VI, ties.method = "min"),
total_rank = rank_sil + rank_dunn + rank_corrected_rand + rank_vi
)
best_overall_bd <- results_ranked_bd %>% arrange(total_rank) %>% head(1)
print(best_overall_bd)## Method k Silhouette Dunn Corrected_Rand VI rank_sil rank_dunn
## 1 ward.D2 2 0.6245096 0.09556396 0.9602 0.1120031 1 11
## rank_corrected_rand rank_vi total_rank
## 1 2 2 16In questo caso la combinazione migliore è data da metodo di linkage ‘ward.D2’ e numero cluster pari a 2. Notiamo che in questo caso che il valore Silhouette 0.6245 è più alto rispetto agli altri due approcci (0.431 per PCA e 0.375 usando tutte le variabili), il che indica che i punti all’interno dei cluster sono ben raggruppati e distinti l’uno dall’altro. Ciò supporta l’idea iniziale che queste due variabili siano in grado di dividere bene il dataset. Il Dunn index risulta più basso rispetto a quello ottenuto con tutte le variabili (0.211) o con la PCA (0.106), il che potrebbe significare che alcuni cluster risultano più vicini o con una certa dispersione interna. Il Corrected Rand 0.9602 rimane molto elevato, segnalando una forte coerenza con la struttura di riferimento attesa. L’indice VI è basso, confermando una buona similarità o una discreta aderenza a una struttura di riferimento.
Visualizziamo come il dataset è stato diviso in cluster.
final_bd <- eclust(banknote_bd_scaled, "hclust", k = 2, hc_metric = "euclidean", hc_method = "ward.D2", graph = FALSE)
fviz_cluster(final_bd, geom = "point", ellipse.type = "norm", palette="jco", ggtheme = theme_minimal())
3.2 Clustering Gerarchico Divisivo
Il clustering divisivo è una tecnica gerarchica “top-down” che parte da un unico cluster contenente l’intero dataset e procede a suddividerlo progressivamente in cluster più piccoli. A differenza del clustering agglomerativo, che inizia con ciascun punto come cluster individuale e li unisce gradualmente, il metodo divisivo si concentra sulla separazione iterativa dei dati in gruppi distinti. Ad ogni suddivisione, si cerca di massimizzare la differenza (o la distanza) tra i cluster risultanti e minimizzare la varianza all’interno di ciascun cluster. Il risultato finale è una struttura ad albero (il dendrogramma) che rappresenta le varie divisioni e il livello di dettaglio raggiunto, offrendo una visione della struttura gerarchica dei dati.
Clustering Gerarchico Divisivo con tutte le variabili
Applichiamo il clusering gerarchico divisivo sul nostro dataset Banknote con tutte le variabili. Utilizziamo l’algoritmo DIANA (DIvisive Analysis Clustering) sul dataset standardizzato (scaling) in precedenza e distanza euclidea per calcolare la distanza tra le osservazioni. La funzione ‘eclust’, se non si specifica il parametro k, calcola automaticamente la statistica del gap per stimare il numero ottimale di cluster. Fornisce inoltre informazioni sulla silhouette per tutti i metodi di partizionamento e per il clustering gerarchico.
diana_clustering <- eclust(df, "diana", stand = FALSE, hc_metric = "euclidean")
fviz_dend(diana_clustering, cex = 0.6, show_labels = FALSE, as.ggplot = TRUE)
Osserviamo che la funzione eclust ha stabilito che il numero ottimale di cluster è 3. Proviamo a calcolare la migliore combinazione di numero di cluster con una leggera modifica del codice visto prima.
evaluate_diana_clustering <- function(df, distanze, status, min_k = 2, max_k = 10) {
results_all <- data.frame( k = numeric(),
Silhouette = numeric(),
Dunn = numeric(),
Corrected_Rand = numeric(),
VI = numeric(),
stringsAsFactors = FALSE)
diana_obj <- diana(distanze)
hclust_obj <- as.hclust(diana_obj)
for(k in min_k:max_k) {
cluster_cut <- cutree(hclust_obj, k = k)
sil <- silhouette(cluster_cut, distanze)
avg_sil <- mean(sil[, 3])
dunn_index <- intCriteria(as.matrix(df), cluster_cut, "Dunn")$dunn
clust_stats <- cluster.stats(distanze, status, cluster_cut)
corrected_rand <- clust_stats$corrected.rand
vi_val <- clust_stats$vi
results_all <- rbind(results_all, data.frame(k = k,
Silhouette = avg_sil,
Dunn = dunn_index,
Corrected_Rand = corrected_rand,
VI = vi_val,
stringsAsFactors = FALSE))
}
return(results_all)
}
results_all_diana <- evaluate_diana_clustering(df, distanze, status)
results_all_diana <- results_all_diana %>%
mutate(
rank_sil = rank(-Silhouette, ties.method = "min"),
rank_dunn = rank(-Dunn, ties.method = "min"),
rank_corrected_rand = rank(-Corrected_Rand, ties.method = "min"),
rank_vi = rank(VI, ties.method = "min"),
total_rank = rank_sil + rank_dunn + rank_corrected_rand + rank_vi
)
best_overall_diana <- results_all_diana %>% arrange(total_rank) %>% head(1)
print(best_overall_diana)## k Silhouette Dunn Corrected_Rand VI rank_sil rank_dunn
## 1 2 0.3723112 0.1819793 0.9406015 0.1540906 1 1
## rank_corrected_rand rank_vi total_rank
## 1 1 1 4Anche in questo caso il risultato suggerisce che il miglior numero di cluster del il dataset è 2. I ranghi di Silhouette, Dunn, Corrected Rand e VI sono tutti pari a 1, il che indica che DIANA con k = 2 ha ottenuto i migliori risultati per queste metriche rispetto agli altri valori di k esplorati. Con un clustering gerarchico divisivo su tutte le variabili numeriche del dataset, con numero di cluster pari a 2 otteniamo: - Silhouette pari a 0.3723. Valore moderato, indica una separazione discreta tra i cluster. - Indice di Dunn pari a 0.18198. Valore piuttosto basso, suggerisce che la separazione tra i cluster non è eccezionale. - Corrected Rand Index pari a 0.9406. Valore alto, indica una buona concordanza con la partizione ideale. - Variational Information (VI) pari a 0.1541. Valore basso, suggerisce una buona corrispondenza con la partizione di riferimento.
Visualizziamo come il dataset è stato diviso in cluster.
final_diana_clustering <- eclust(df, "diana", stand = FALSE, k = 2, hc_metric = "euclidean")
fviz_cluster(final_diana_clustering, geom = "point", ellipse.type = "norm", palette="jco", ggtheme = theme_minimal())
Clustering Gerarchico Divisivo con le componenti principali
Applichiamo il clusering gerarchico divisivo sul nostro dataset considerando solo le 3 componenti principali.
results_pca_diana <- evaluate_diana_clustering(banknote_pca_num, dist(banknote_pca_num), status)
results_pca_diana <- results_pca_diana %>%
mutate(
rank_sil = rank(-Silhouette, ties.method = "min"),
rank_dunn = rank(-Dunn, ties.method = "min"),
rank_corrected_rand = rank(-Corrected_Rand, ties.method = "min"),
rank_vi = rank(VI, ties.method = "min"),
total_rank = rank_sil + rank_dunn + rank_corrected_rand + rank_vi
)
best_pca_diana <- results_pca_diana %>% arrange(total_rank) %>% head(1)
print(best_pca_diana)## k Silhouette Dunn Corrected_Rand VI rank_sil rank_dunn
## 1 2 0.4393924 0.09331653 0.8456292 0.2819728 1 1
## rank_corrected_rand rank_vi total_rank
## 1 1 1 4Ancora una volta viene indicato 2 come numero ideali di cluster.
I risultati ottenuti applicando PCA prima di eseguire il clustering DIANA per k = 2 sono i seguenti:
Silhouette pari a 0.4394, leggermente migliore rispetto a quando non è stata applicata la PCA, suggerendo una separazione più chiara tra i cluster.
Indice di Dunn pari a 0.0933, piuttosto basso, indicante una separazione limitata tra i cluster.
Corrected Rand Index pari a 0.8456, buona concordanza con la partizione ideale, ma inferiore a quando non si applica la PCA.
Variational Information (VI) pari a 0.2820, un po’ più alta rispetto al caso senza PCA, indicando una leggera differenza rispetto alla partizione di riferimento.
L’applicazione della PCA prima del clustering DIANA porta a un miglioramento nella separazione dei cluster, come evidenziato dal valore di Silhouette più alto (0.4394 rispetto a 0.3723), ma l’indice di Dunn rimane basso, suggerendo che i cluster non sono completamente separati. Il Corrected Rand Index e il Variational Information (VI) mostrano una discreta corrispondenza con una partizione ideale, ma comunque inferiore rispetto ai risultati ottenuti senza PCA.
Visualizziamo come il dataset è stato diviso in cluster.
final_diana_clustering_pca <- eclust(banknote_pca_num, "diana", stand = FALSE, k = 2, hc_metric = "euclidean")
clusters_diana <- final_diana_clustering_pca$cluster
scatterplot3d(banknote_pca_num,
color = clusters_diana,
pch = 19,
main = "DIANA Clustering su PCA (PC1, PC2, PC3)",
xlab = "PC1", ylab = "PC2", zlab = "PC3")
legend("topright", legend = unique(clusters_diana), col = unique(clusters_diana), pch = 19, title = "Cluster")
Clustering Gerarchico Divisivo con dataset ridotto
Applichiamo il clusering gerarchico divisivo sul nostro dataset ridotto alle sole due variabili ‘Bottom’ e ‘Diagonal’.
results_bd_diana <- evaluate_diana_clustering(banknote_bd_scaled, dist(banknote_bd_scaled), status)
results_bd_diana <- results_bd_diana %>%
mutate(
rank_sil = rank(-Silhouette, ties.method = "min"),
rank_dunn = rank(-Dunn, ties.method = "min"),
rank_corrected_rand = rank(-Corrected_Rand, ties.method = "min"),
rank_vi = rank(VI, ties.method = "min"),
total_rank = rank_sil + rank_dunn + rank_corrected_rand + rank_vi
)
best_bd_diana <- results_bd_diana %>% arrange(total_rank) %>% head(1)
print(best_bd_diana)## k Silhouette Dunn Corrected_Rand VI rank_sil rank_dunn
## 1 2 0.6241219 0.07897225 0.9406015 0.1540906 2 2
## rank_corrected_rand rank_vi total_rank
## 1 1 1 6I risultati ottenuti applicando DIANA al dataset composto solo dalle variabili ‘Bottom’ e ‘Diagonal’ per k = 2 sono i seguenti:
Silhouette pari a 0.6241, valore buono, indica una separazione più marcata tra i cluster rispetto ai precedenti risultati.
Indice di Dunn pari a 0.07897, valore ancora piuttosto basso, suggerisce che la separazione tra i cluster non è ottimale, ma non peggiora rispetto ai casi precedenti.
Corrected Rand Index pari a 0.9406, valore molto alto, indica che il clustering ottenuto con DIANA è molto simile alla partizione ideale.
Variational Information (VI) pari a 0.1541, valore basso, indica una buona corrispondenza con la partizione di riferimento.
Applicare DIANA solo sulle variabili bottom e diagonal migliora il risultato della separazione dei cluster, come evidenziato dal valore più alto di Silhouette (0.6241). Tuttavia, l’indice di Dunn rimane basso, suggerendo che la separazione tra i cluster non è ancora ideale. La concordanza con la partizione ideale (Corrected Rand Index) è molto alta, il che indica che il clustering rispecchia bene una divisione naturale del dataset.
Visualizziamo come il dataset è stato diviso in cluster.
final_diana_clustering_bd <- eclust(banknote_bd_scaled, "diana", stand = FALSE, k = 2, hc_metric = "euclidean")
fviz_cluster(final_diana_clustering_bd, geom = "point", ellipse.type = "norm", palette="jco", ggtheme = theme_minimal())
Vantaggi e Svantaggi del Clustering Gerarchico
TODO #### Vantaggi:
Non richiede la specifica del numero di cluster a priori.
Produce una rappresentazione gerarchica che può essere utile per comprendere le relazioni tra i dati.
Adatto a dataset di piccole e medie dimensioni.Svantaggi:
TODO
Sensibile alla scelta della misura di distanza e del metodo di linkage.
Può essere computazionalmente costoso per dataset di grandi dimensioni. (quello divisivo)
Una volta effettuata una fusione, non è possibile separare i cluster, rendendo l'algoritmo meno flessibile in caso di errori.3.2 Clustering Partizionale: K-means e K-medoids (PAM)
3.2.1 K-Means
3.2.1.1 Esecuzione dell’algoritmo con tutte le variabili
In questo capitolo i cluster vengono creati utilizzando l’algoritmo k-means, dopo aver scalato le variabili numeriche del dataset.
library(factoextra)
df <- scale(banknote[, -1])Successivamente si determina il numero ottimale di cluster da utilizzare nell’algoritmo, impiegando il metodo Elbow.
fviz_nbclust(df, kmeans, method = "wss") +
geom_vline(xintercept = 4, linetype = 2)+
labs(subtitle = "Elbow method")
Il metodo Silhouette.
fviz_nbclust(df, kmeans, method = "silhouette")+
labs(subtitle = "Silhouette method")
Ed infine la GAP analysis:
set.seed(123)
fviz_nbclust(df, kmeans, nstart = 25, method = "gap_stat", nboot = 50)+
labs(subtitle = "Gap statistic method")
Il primo metodo indica quattro come numero ottimale di cluster, il secondo ne suggerisce due mentre la gap statistic tre. Poiché il dataset è composto da 200 osservazioni, equamente suddivise tra banconote autentiche e contraffatte, si sceglie di impostare il numero di cluster a due e si procede con l’esecuzione dell’algoritmo.
Per prima cosa si osserva una matrice di grafici a dispersione, utile per analizzare le relazioni tra le variabili numeriche del dataset.
set.seed(123)
km.res <- kmeans(df, 2, nstart = 25)
print(km.res)## K-means clustering with 2 clusters of sizes 108, 92
##
## Cluster means:
## Length Left Right Bottom Top Diagonal
## 1 -0.1221604 0.5482838 0.6163644 0.6788357 0.5494346 -0.8005001
## 2 0.1434057 -0.6436375 -0.7235582 -0.7968941 -0.6449884 0.9397176
##
## Clustering vector:
## [1] 1 2 2 2 2 1 2 2 2 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [38] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 1 2 2 2 2
## [75] 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 1 1 1 1 1 1 1 1 1 1 1
## [112] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [149] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [186] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##
## Within cluster sum of squares by cluster:
## [1] 396.3549 304.8505
## (between_SS / total_SS = 41.3 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss" "tot.withinss"
## [6] "betweenss" "size" "iter" "ifault"aggregate(banknote, by=list(cluster=km.res$cluster), mean)## cluster Status Length Left Right Bottom Top Diagonal
## 1 1 NA 214.85 130.3194 130.2056 10.398148 11.09167 139.5611
## 2 2 NA 214.95 129.8891 129.6641 8.266304 10.13261 141.5663dd <- cbind(banknote, cluster = km.res$cluster)
head(dd)## Status Length Left Right Bottom Top Diagonal cluster
## 1 genuine 214.8 131.0 131.1 9.0 9.7 141.0 1
## 2 genuine 214.6 129.7 129.7 8.1 9.5 141.7 2
## 3 genuine 214.8 129.7 129.7 8.7 9.6 142.2 2
## 4 genuine 214.8 129.7 129.6 7.5 10.4 142.0 2
## 5 genuine 215.0 129.6 129.7 10.4 7.7 141.8 2
## 6 genuine 215.7 130.8 130.5 9.0 10.1 141.4 1km.res$cluster## [1] 1 2 2 2 2 1 2 2 2 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [38] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 1 2 2 2 2
## [75] 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 1 1 1 1 1 1 1 1 1 1 1
## [112] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [149] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [186] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1km.res$size## [1] 108 92km.res$centers## Length Left Right Bottom Top Diagonal
## 1 -0.1221604 0.5482838 0.6163644 0.6788357 0.5494346 -0.8005001
## 2 0.1434057 -0.6436375 -0.7235582 -0.7968941 -0.6449884 0.9397176cl <- km.res$cluster
pairs(df, gap=0, pch=cl, col=c("red", "blue")[cl])
TODO: commentare
Ora si analizza la distribuzione dei due cluster.
fviz_cluster(km.res,
data = df,
palette = c("red", "blue"),
ellipse.type = "euclid",
star.plot = TRUE,
repel = TRUE,
ggtheme = theme_minimal()
)
TODO: aggiungere eclust() (che utilizzare gap statistics per calcolare k)
3.2.1.2 Esecuzione dell’algoritmo con solamente le variabili bottom e diagonal
In questo capitolo i cluster vengono creati utilizzando l’algoritmo k-means, dopo aver scalato le variabili numeriche del dataset.
df_ridotto <- df[, c("Bottom", "Diagonal")]TODO: Successivamente si determina il numero ottimale di cluster da utilizzare nell’algoritmo, impiegando il metodo Elbow.
fviz_nbclust(df_ridotto, kmeans, method = "wss") +
geom_vline(xintercept = 4, linetype = 2)+
labs(subtitle = "Elbow method")
TODO: Il metodo Silhouette.
fviz_nbclust(df_ridotto, kmeans, method = "silhouette")+
labs(subtitle = "Silhouette method")
TODO: Ed infine la GAP analysis:
set.seed(123)
fviz_nbclust(df_ridotto, kmeans, nstart = 25, method = "gap_stat", nboot = 500)+
labs(subtitle = "Gap statistic method")
TODO: Il primo metodo indica quattro come numero ottimale di cluster, il secondo ne suggerisce due mentre la gap statistic tre. Poiché il dataset è composto da 200 osservazioni, equamente suddivise tra banconote autentiche e contraffatte, si sceglie di impostare il numero di cluster a due e si procede con l’esecuzione dell’algoritmo.
TODO: Per prima cosa si osserva una matrice di grafici a dispersione, utile per analizzare le relazioni tra le variabili numeriche del dataset.
set.seed(123)
km.res <- kmeans(df_ridotto, 2, nstart = 25)
print(km.res)## K-means clustering with 2 clusters of sizes 99, 101
##
## Cluster means:
## Bottom Diagonal
## 1 0.7829035 -0.9035032
## 2 -0.7674005 0.8856120
##
## Clustering vector:
## [1] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [38] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [75] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 2 1 1 1 1 1 1 1 1
## [112] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [149] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [186] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##
## Within cluster sum of squares by cluster:
## [1] 81.94094 35.86830
## (between_SS / total_SS = 70.4 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss" "tot.withinss"
## [6] "betweenss" "size" "iter" "ifault"aggregate(banknote, by=list(cluster=km.res$cluster), mean)## cluster Status Length Left Right Bottom Top Diagonal
## 1 1 NA 214.8222 130.3000 130.1939 10.548485 11.12727 139.4424
## 2 2 NA 214.9683 129.9465 129.7238 8.308911 10.18317 141.5040dd <- cbind(banknote, cluster = km.res$cluster)
head(dd)## Status Length Left Right Bottom Top Diagonal cluster
## 1 genuine 214.8 131.0 131.1 9.0 9.7 141.0 2
## 2 genuine 214.6 129.7 129.7 8.1 9.5 141.7 2
## 3 genuine 214.8 129.7 129.7 8.7 9.6 142.2 2
## 4 genuine 214.8 129.7 129.6 7.5 10.4 142.0 2
## 5 genuine 215.0 129.6 129.7 10.4 7.7 141.8 2
## 6 genuine 215.7 130.8 130.5 9.0 10.1 141.4 2km.res$cluster## [1] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [38] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [75] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 2 1 1 1 1 1 1 1 1
## [112] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [149] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [186] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1km.res$size## [1] 99 101km.res$centers## Bottom Diagonal
## 1 0.7829035 -0.9035032
## 2 -0.7674005 0.8856120cl <- km.res$cluster
pairs(df_ridotto, gap=0, pch=cl, col=c("red", "blue")[cl])
TODO: commentare
TODO: Ora si analizza la distribuzione dei due cluster.
fviz_cluster(km.res,
data = df_ridotto,
palette = c("red", "blue"),
ellipse.type = "euclid",
star.plot = TRUE,
repel = TRUE,
ggtheme = theme_minimal()
)
3.2.1.3 Esecuzione dell’algoritmo con solamente le variabili scelte tramite la PCA
In questo capitolo i cluster vengono creati utilizzando l’algoritmo k-means, dopo aver scalato le variabili numeriche del dataset.
df_pca <- scale(banknote_pca[, -1])TODO: Successivamente si determina il numero ottimale di cluster da utilizzare nell’algoritmo, impiegando il metodo Elbow.
fviz_nbclust(df_pca, kmeans, method = "wss") +
geom_vline(xintercept = 4, linetype = 2)+
labs(subtitle = "Elbow method")
TODO: Il metodo Silhouette.
fviz_nbclust(df_pca, kmeans, method = "silhouette")+
labs(subtitle = "Silhouette method")
TODO: Ed infine la GAP analysis:
set.seed(123)
fviz_nbclust(df_pca, kmeans, nstart = 25, method = "gap_stat", nboot = 500)+
labs(subtitle = "Gap statistic method")
TODO: Il primo metodo indica quattro come numero ottimale di cluster, il secondo ne suggerisce due mentre la gap statistic tre. Poiché il dataset è composto da 200 osservazioni, equamente suddivise tra banconote autentiche e contraffatte, si sceglie di impostare il numero di cluster a due e si procede con l’esecuzione dell’algoritmo.
TODO: Per prima cosa si osserva una matrice di grafici a dispersione, utile per analizzare le relazioni tra le variabili numeriche del dataset.
set.seed(123)
km.res <- kmeans(df_pca, 2, nstart = 25)
print(km.res)## K-means clustering with 2 clusters of sizes 100, 100
##
## Cluster means:
## PCA1 PCA2 PCA3
## 1 -0.8574889 0.3668574 -0.02482513
## 2 0.8574889 -0.3668574 0.02482513
##
## Clustering vector:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
## 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80
## 1 1 1 1 1 1 1 1 1 2 2 1 1 1 1 1 1 1 1 1
## 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120
## 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2
## 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140
## 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160
## 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180
## 1 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2
## 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200
## 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
##
## Within cluster sum of squares by cluster:
## [1] 221.1989 201.7035
## (between_SS / total_SS = 29.2 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss" "tot.withinss"
## [6] "betweenss" "size" "iter" "ifault"aggregate(banknote, by=list(cluster=km.res$cluster), mean)## cluster Status Length Left Right Bottom Top Diagonal
## 1 1 NA 215.008 129.945 129.713 8.307 10.183 141.466
## 2 2 NA 214.784 130.298 130.200 10.528 11.118 139.501dd <- cbind(banknote, cluster = km.res$cluster)
head(dd)## Status Length Left Right Bottom Top Diagonal cluster
## 1 genuine 214.8 131.0 131.1 9.0 9.7 141.0 2
## 2 genuine 214.6 129.7 129.7 8.1 9.5 141.7 1
## 3 genuine 214.8 129.7 129.7 8.7 9.6 142.2 1
## 4 genuine 214.8 129.7 129.6 7.5 10.4 142.0 1
## 5 genuine 215.0 129.6 129.7 10.4 7.7 141.8 1
## 6 genuine 215.7 130.8 130.5 9.0 10.1 141.4 1km.res$cluster## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
## 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80
## 1 1 1 1 1 1 1 1 1 2 2 1 1 1 1 1 1 1 1 1
## 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120
## 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2
## 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140
## 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160
## 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180
## 1 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2
## 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200
## 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2km.res$size## [1] 100 100km.res$centers## PCA1 PCA2 PCA3
## 1 -0.8574889 0.3668574 -0.02482513
## 2 0.8574889 -0.3668574 0.02482513cl <- km.res$cluster
pairs(df_ridotto, gap=0, pch=cl, col=c("red", "blue")[cl])
TODO: commentare
TODO: Ora si analizza la distribuzione dei due cluster.
fviz_cluster(km.res,
data = df_pca,
palette = c("red", "blue"),
ellipse.type = "euclid",
star.plot = TRUE,
repel = TRUE,
ggtheme = theme_minimal()
)
3.2.2 K-Medoids
3.2.2.1 Esecuzione dell’algoritmo sul dataset selezionando tutte le variabili
In questo capitolo i cluster vengono creati utilizzando l’algoritmo k-medoids con k impostato a 2.
Come in precedenza, si analizzano prima i grafici a dispersione e successivamente i due cluster.
library(cluster)
head(df, n = 3)## Length Left Right Bottom Top Diagonal
## [1,] -0.2549435 2.433346 2.829942 -0.2890067 -1.183765 0.4482473
## [2,] -0.7860757 -1.167507 -0.634788 -0.9120152 -1.432847 1.0557460
## [3,] -0.2549435 -1.167507 -0.634788 -0.4966762 -1.308306 1.4896737pam.res <- pam(df, 2)
print(pam.res)## Medoids:
## ID Length Left Right Bottom Top Diagonal
## [1,] 193 -0.5205096 0.4944248 0.6026155 0.9570103 0.5598130 -1.113892
## [2,] 47 -0.2549435 -0.6135300 -0.6347880 -0.7735689 -0.5610584 0.882175
## Clustering vector:
## [1] 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [38] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2
## [75] 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1
## [112] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [149] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [186] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## Objective function:
## build swap
## 1.842553 1.788867
##
## Available components:
## [1] "medoids" "id.med" "clustering" "objective" "isolation"
## [6] "clusinfo" "silinfo" "diss" "call" "data"dd <- cbind(banknote, cluster = pam.res$cluster)
head(dd, n = 8)## Status Length Left Right Bottom Top Diagonal cluster
## 1 genuine 214.8 131.0 131.1 9.0 9.7 141.0 1
## 2 genuine 214.6 129.7 129.7 8.1 9.5 141.7 2
## 3 genuine 214.8 129.7 129.7 8.7 9.6 142.2 2
## 4 genuine 214.8 129.7 129.6 7.5 10.4 142.0 2
## 5 genuine 215.0 129.6 129.7 10.4 7.7 141.8 2
## 6 genuine 215.7 130.8 130.5 9.0 10.1 141.4 2
## 7 genuine 215.5 129.5 129.7 7.9 9.6 141.6 2
## 8 genuine 214.5 129.6 129.2 7.2 10.7 141.7 2pam.res$medoids## Length Left Right Bottom Top Diagonal
## [1,] -0.5205096 0.4944248 0.6026155 0.9570103 0.5598130 -1.113892
## [2,] -0.2549435 -0.6135300 -0.6347880 -0.7735689 -0.5610584 0.882175head(pam.res$clustering)## [1] 1 2 2 2 2 2cl <- pam.res$clustering
pairs(df, gap=0, pch=cl, col=c("#00AFBB", "#FC4E07")[cl])
Come previsto, i cluster risultano molto simili.
fviz_cluster(pam.res,
palette = c("#00AFBB", "#FC4E07"),
ellipse.type = "t",
repel = TRUE,
ggtheme = theme_classic()
)
3.2.2.2 Esecuzione dell’algoritmo con solamente le variabili bottom e diagonal
df_ridotto <- df[, c("Bottom", "Diagonal")]
pam.res <- pam(df_ridotto, 2)
print(pam.res)## Medoids:
## ID Bottom Diagonal
## [1,] 47 -0.7735689 0.8821750
## [2,] 185 0.8877871 -0.8535357
## Clustering vector:
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [38] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [75] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 1 2 2 2 2 2 2 2 2
## [112] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [149] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## [186] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## Objective function:
## build swap
## 0.8831314 0.6414212
##
## Available components:
## [1] "medoids" "id.med" "clustering" "objective" "isolation"
## [6] "clusinfo" "silinfo" "diss" "call" "data"dd <- cbind(banknote, cluster = pam.res$cluster)
head(dd, n = 8)## Status Length Left Right Bottom Top Diagonal cluster
## 1 genuine 214.8 131.0 131.1 9.0 9.7 141.0 1
## 2 genuine 214.6 129.7 129.7 8.1 9.5 141.7 1
## 3 genuine 214.8 129.7 129.7 8.7 9.6 142.2 1
## 4 genuine 214.8 129.7 129.6 7.5 10.4 142.0 1
## 5 genuine 215.0 129.6 129.7 10.4 7.7 141.8 1
## 6 genuine 215.7 130.8 130.5 9.0 10.1 141.4 1
## 7 genuine 215.5 129.5 129.7 7.9 9.6 141.6 1
## 8 genuine 214.5 129.6 129.2 7.2 10.7 141.7 1pam.res$medoids## Bottom Diagonal
## [1,] -0.7735689 0.8821750
## [2,] 0.8877871 -0.8535357head(pam.res$clustering)## [1] 1 1 1 1 1 1cl <- pam.res$clustering
pairs(df, gap=0, pch=cl, col=c("#00AFBB", "#FC4E07")[cl])
Come previsto, i cluster risultano molto simili.
fviz_cluster(pam.res,
palette = c("#00AFBB", "#FC4E07"),
ellipse.type = "t",
repel = TRUE,
ggtheme = theme_classic()
)
3.2.2.3 Esecuzione dell’algoritmo con solamente le variabili scelte tramite la PCA
df_pca <- scale(banknote_pca[, -1])
pam.res <- pam(df_pca, 2)
print(pam.res)## Medoids:
## ID PCA1 PCA2 PCA3
## 198 198 0.8816755 -0.1289878 0.19323058
## 21 21 -1.1983358 0.0749874 0.05222386
## Clustering vector:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
## 1 2 2 2 2 1 2 2 2 1 1 2 2 2 2 2 2 2 2 2
## 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
## 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
## 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
## 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2
## 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80
## 2 2 2 2 2 1 2 2 2 1 2 2 2 2 2 2 2 2 2 2
## 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
## 2 2 2 2 1 2 2 2 2 2 2 2 2 2 2 2 1 2 2 2
## 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200
## 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## Objective function:
## build swap
## 1.370309 1.329691
##
## Available components:
## [1] "medoids" "id.med" "clustering" "objective" "isolation"
## [6] "clusinfo" "silinfo" "diss" "call" "data"dd <- cbind(banknote, cluster = pam.res$cluster)
head(dd, n = 8)## Status Length Left Right Bottom Top Diagonal cluster
## 1 genuine 214.8 131.0 131.1 9.0 9.7 141.0 1
## 2 genuine 214.6 129.7 129.7 8.1 9.5 141.7 2
## 3 genuine 214.8 129.7 129.7 8.7 9.6 142.2 2
## 4 genuine 214.8 129.7 129.6 7.5 10.4 142.0 2
## 5 genuine 215.0 129.6 129.7 10.4 7.7 141.8 2
## 6 genuine 215.7 130.8 130.5 9.0 10.1 141.4 1
## 7 genuine 215.5 129.5 129.7 7.9 9.6 141.6 2
## 8 genuine 214.5 129.6 129.2 7.2 10.7 141.7 2pam.res$medoids## PCA1 PCA2 PCA3
## 198 0.8816755 -0.1289878 0.19323058
## 21 -1.1983358 0.0749874 0.05222386head(pam.res$clustering)## 1 2 3 4 5 6
## 1 2 2 2 2 1cl <- pam.res$clustering
pairs(df, gap=0, pch=cl, col=c("#00AFBB", "#FC4E07")[cl])
Come previsto, i cluster risultano molto simili.
fviz_cluster(pam.res,
palette = c("#00AFBB", "#FC4E07"),
ellipse.type = "t",
repel = TRUE,
ggtheme = theme_classic()
)
3.1.4 Misture di Gaussiane
In questa sezione introduciamo l’applicazione del clustering basato
su misture di gaussiane, utilizzando la funzione Mclust del
pacchetto mclust. Questo approccio assume che i dati
siano generati da una combinazione di distribuzioni gaussiane e, tramite
l’algoritmo EM (Expectation-Maximization), stima i parametri di ciascuna
componente. Il vantaggio principale di questo metodo è la capacità di
modellare cluster con forme ellissoidali, permettendo di gestire anche
cluster con sovrapposizioni. Inoltre, Mclust effettua una
selezione automatica del numero ottimale di cluster e del modello di
covarianza (parsimonioso) in base al criterio BIC, fornendo
un’indicazione della bontà dell’adattamento e della complessità del
modello.
library(mclust)
data(banknote)
dati <- banknote[, -1]
dati_scaled <- scale(dati)
gmm_result <- Mclust(dati_scaled)
summary(gmm_result)## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm
## ----------------------------------------------------
##
## Mclust VEE (ellipsoidal, equal shape and orientation) model with 4 components:
##
## log-likelihood n df BIC ICL
## -1191.595 200 51 -2653.405 -2666.898
##
## Clustering table:
## 1 2 3 4
## 16 99 47 38plot(gmm_result, what = "BIC")
plot(gmm_result, what = "classification")
plot(gmm_result, what = "uncertainty")
plot(gmm_result, what = "density")
Nonostante sappiamo che il dataset banknote presenta due classi reali
(ad esempio, banconote fraudolente e non fraudolente), l’analisi con
Mclust ha individuato come miglior modello un numero di cluster pari a
3. Inoltre, il modello selezionato è di tipo VEE (ovvero ellissoidale,
con forma variabile, ma con volume ed orientamento uguali), in base al
valore più alto di BIC evidenziato nel plot. È interessante notare come
il BIC per 2 cluster risulti significativamente inferiore rispetto a
quello per un numero maggiore di cluster, probabilmente perché sono
state utilizzate tutte le variabili del dataset per costruire i
cluster.
Confrontiamo esplicitamente i valori di BIC e ICL per soluzioni con 2 e
3 cluster:
library(mclust)
data(banknote)
dati <- banknote[, -1]
dati_scaled <- scale(dati)
gmm_result2 <- Mclust(dati_scaled, G=2)
summary(gmm_result2)## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm
## ----------------------------------------------------
##
## Mclust EEV (ellipsoidal, equal volume and shape) model with 2 components:
##
## log-likelihood n df BIC ICL
## -1260.823 200 49 -2781.264 -2781.417
##
## Clustering table:
## 1 2
## 101 99gmm_result3 <- Mclust(dati_scaled, G=3)
summary(gmm_result3)## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm
## ----------------------------------------------------
##
## Mclust VEE (ellipsoidal, equal shape and orientation) model with 3 components:
##
## log-likelihood n df BIC ICL
## -1212.796 200 43 -2653.42 -2653.74
##
## Clustering table:
## 1 2 3
## 16 99 85Dai sommari ottenuti si osserva che entrambi i modelli, per k = 2 e k
= 3, forniscono valori che suggeriscono una soluzione migliore per 3
cluster. I valori del BIC e, conseguentemente, anche quelli dell’ICL
(che include una penalizzazione per l’incertezza nelle assegnazioni),
indicano che la soluzione a 3 cluster è preferibile rispetto a quella a
2 cluster.
Valutiamo l’effetto della scelta delle variabili sul modello.
Utilizziamo invece solo le due variabili “Bottom” e “Diagonal”,
individuate in precedenza dagli scatterplot come le più
discriminanti:
library(mclust)
data(banknote)
dati <- banknote[, c("Bottom", "Diagonal")]
dati_scaled <- scale(dati)
gmm_result <- Mclust(dati_scaled)
summary(gmm_result)## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm
## ----------------------------------------------------
##
## Mclust EVE (ellipsoidal, equal volume and orientation) model with 3 components:
##
## log-likelihood n df BIC ICL
## -387.8563 200 13 -844.5908 -849.2887
##
## Clustering table:
## 1 2 3
## 100 16 84plot(gmm_result, what = "BIC")
plot(gmm_result, what = "classification")
plot(gmm_result, what = "uncertainty")
plot(gmm_result, what = "density")
Anche in questo caso si nota che il valore di BIC è notevolmente più
alto (migliore) rispetto al caso in cui sono state usate tutte le
variabili. L’algoritmo suggerisce ancora un numero ottimale di cluster
pari a 3, ma questa volta ha selezionato come modello parsimonioso il
modello EVE (ellissoidale, con volume uguale e orientamento variabile).
È importante osservare che in uno dei cluster sono presenti solamente 16
osservazioni, il che potrebbe indicare un gruppo di dati che,
probabilmente, sono state erroneamente assegnate a un cluster
separato.
Verifichiamo ora se i valori di ICL sono in accordo con quelli del BIC
confrontando esplicitamente le soluzioni per 2 e 3 cluster:
library(mclust)
data(banknote)
dati <- banknote[, c("Bottom", "Diagonal")]
dati_scaled <- scale(dati)
gmm_result2 <- Mclust(dati_scaled, G=2)
summary(gmm_result2)## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm
## ----------------------------------------------------
##
## Mclust VVE (ellipsoidal, equal orientation) model with 2 components:
##
## log-likelihood n df BIC ICL
## -417.7539 200 10 -888.491 -890.1796
##
## Clustering table:
## 1 2
## 99 101gmm_result3 <- Mclust(dati_scaled, G=3)
summary(gmm_result3)## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm
## ----------------------------------------------------
##
## Mclust EVE (ellipsoidal, equal volume and orientation) model with 3 components:
##
## log-likelihood n df BIC ICL
## -387.8563 200 13 -844.5908 -849.2887
##
## Clustering table:
## 1 2 3
## 100 16 84In questo caso, i valori di BIC e ICL sono concordi: per un numero di
cluster pari a 3 si ottiene un BIC di circa -845 e un ICL di -849,
mentre per 2 cluster il BIC risulta attorno a -888 e l’ICL a -890.
Questo dimostra che anche l’ICL, che penalizza ulteriormente le
assegnazioni incerte, suggerisce che il modello ottimale è quello con 3
cluster.
Adesso, procederemo ad utilizzare solamente le variabili estratte
tramite la tecnica PCA.
dati <- banknote_pca[, -1]
dati_scaled <- scale(dati)
gmm_result <- Mclust(dati_scaled)
summary(gmm_result)## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm
## ----------------------------------------------------
##
## Mclust EEV (ellipsoidal, equal volume and shape) model with 2 components:
##
## log-likelihood n df BIC ICL
## -764.5395 200 16 -1613.852 -1616.362
##
## Clustering table:
## 1 2
## 102 98plot(gmm_result, what = "BIC")
plot(gmm_result, what = "classification")
plot(gmm_result, what = "uncertainty")
plot(gmm_result, what = "density")
In questo caso si nota che il valore di BIC è più alto (migliore)
rispetto al caso in cui sono state usate tutte le variabili ma è più
basso del caso in cui sono state usate solamente le due variabili
selezionate. Questa volta l’algoritmo suggerisce un numero ottimale di
cluster pari a 2, selezionando come modello parsimonioso il modello EEV
(ellissoidale, con volume e orientamento uguale).
Verifichiamo ora se i valori di ICL sono in accordo con quelli del BIC
confrontando esplicitamente le soluzioni per 2 e 3 cluster:
gmm_result2 <- Mclust(dati_scaled, G=2)
summary(gmm_result2)## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm
## ----------------------------------------------------
##
## Mclust EEV (ellipsoidal, equal volume and shape) model with 2 components:
##
## log-likelihood n df BIC ICL
## -764.5395 200 16 -1613.852 -1616.362
##
## Clustering table:
## 1 2
## 102 98gmm_result3 <- Mclust(dati_scaled, G=3)
summary(gmm_result3)## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm
## ----------------------------------------------------
##
## Mclust VEE (ellipsoidal, equal shape and orientation) model with 3 components:
##
## log-likelihood n df BIC ICL
## -764.778 200 19 -1630.224 -1657.822
##
## Clustering table:
## 1 2 3
## 64 98 38Anche in questo caso, i valori di BIC e ICL sono concordi: per un
numero di cluster pari a 3 si ottiene un BIC di circa -1630 e un ICL di
-1657, mentre per 2 cluster il BIC risulta attorno a -1614 e l’ICL a
-1616. Questo dimostra che anche l’ICL suggerisce che il modello
ottimale è quello con 2 cluster.
## 3.1.5. Validazione e visualizzazione dei risultati
TODO: inserire commento introduttivo
3.1.5.3. Validazione con la classe reale
In questa sezione analizzeremo la qualità dei diversi metodi di
clustering applicati al dataset banknote, fissando il numero di cluster
pari al numero di classi reali, ovvero 2. Questo approccio ci permette
di utilizzare tecniche di validazione esterne, dato che disponiamo delle
etichette reali delle osservazioni (variabile ‘Status’).
Per valutare l’efficacia dei cluster identificati, utilizzeremo diverse
metriche, tra cui:
- Adjusted Rand Index (ARI) e Meila Variation Index, che misurano il
grado di concordanza tra il clustering ottenuto e la classificazione
reale. - Metriche di valutazione proprie del machine learning
supervisionato, come accuratezza, specificità, recall, matrice di
confusione e AUC, per interpretare il clustering come un problema di
classificazione binaria.
Verranno valutate le misure di validazione al variare del numero di
variabili utilizzate nel clustering. In particolare, analizzeremo tre
configurazioni differenti:
- Dataset completo: utilizzo di tutte le variabili
disponibili. - Dataset ridotto: utilizzo esclusivo
delle variabili ‘Diagonal’ e ‘Bottom’, identificate come le più
discriminative. - Dataset trasformato: utilizzo delle
variabili ottenute tramite PCA.
3.1.5.3.1. Dataset completo
Nel seguente codice, applichiamo i diversi algoritmi di clustering al dataset originale con i migliori parametri calcolati in precedenza (fissando il numero di cluster a 2) e aggiungiamo al dataframe le colonne corrispondenti alle assegnazioni di ciascun metodo.
library(cluster)
data("banknote")
data <- banknote[, -1]
data_scaled <- scale(data)
set.seed(123)
km_model <- kmeans(data_scaled, centers = 2, nstart = 25)
pam_model <- pam(data_scaled, k = 2, nstart = 25)
hc_model <- hclust(dist(data_scaled, method = "euclidean"), method = "ward.D2")
diana_model <- diana(data_scaled)
gmm_model <- Mclust(data_scaled, G = 2)
banknote$kmeans <- km_model$cluster
banknote$kmedoids <- pam_model$clustering
banknote$hierarchical_agglomerative <- cutree(hc_model, k = 2)
banknote$hierarchical_divisive <- cutree(as.hclust(diana_model), k = 2)
banknote$gmm <- gmm_model$classification
head(banknote)## Status Length Left Right Bottom Top Diagonal kmeans kmedoids
## 1 genuine 214.8 131.0 131.1 9.0 9.7 141.0 1 1
## 2 genuine 214.6 129.7 129.7 8.1 9.5 141.7 2 2
## 3 genuine 214.8 129.7 129.7 8.7 9.6 142.2 2 2
## 4 genuine 214.8 129.7 129.6 7.5 10.4 142.0 2 2
## 5 genuine 215.0 129.6 129.7 10.4 7.7 141.8 2 2
## 6 genuine 215.7 130.8 130.5 9.0 10.1 141.4 1 2
## hierarchical_agglomerative hierarchical_divisive gmm
## 1 1 1 2
## 2 2 2 2
## 3 2 2 2
## 4 2 2 2
## 5 2 2 2
## 6 2 2 2Il codice successivo implementa una funzione per riallineare le assegnazioni dei cluster alle classi reali in modo da rendere i confronti più accurati.
library(caret)
library(pROC)
library(mclust)
bestMapping <- function(true, pred) {
banknote <- banknote %>%
mutate(Status_numeric = recode(Status,
"counterfeit" = 1,
"genuine" = 2))
true <- banknote$Status_numeric
true <- as.character(true)
pred <- as.character(pred)
levels_true <- sort(unique(true))
cm1 <- confusionMatrix(factor(pred, levels = levels_true), factor(true, levels = levels_true))
acc1 <- cm1$overall["Accuracy"]
pred_swap <- ifelse(pred == levels_true[1], levels_true[2], levels_true[1])
cm2 <- confusionMatrix(factor(pred_swap, levels = levels_true), factor(true, levels = levels_true))
acc2 <- cm2$overall["Accuracy"]
if(acc2 > acc1) {
return(list(mapped = pred_swap, cm = cm2))
} else {
return(list(mapped = pred, cm = cm1))
}
}Nel blocco successivo, calcoliamo le metriche di valutazione per ciascun algoritmo di clustering. È importante notare che il valore di AUC ha un significato più chiaro per l’algoritmo GMM (Gaussian Mixture Model), poiché questo metodo fornisce direttamente le probabilità a posteriori di appartenenza ai cluster. Per gli altri algoritmi, che operano in modo hard (assegnando ciascun punto a un solo cluster senza una misura di incertezza), il calcolo dell’AUC è stato forzato attraverso:
- Per K-means e K-medoids, è stato definito uno score basato sulla
distanza dei punti dal centroide o dal medoide del rispettivo
cluster.
- Per gli algoritmi gerarchici, non è stato possibile adottare un criterio analogo per costruire uno score continuo, motivo per cui l’AUC non è stato calcolato (indicata come NA nella tabella).
# K-means
ari_kmeans <- adjustedRandIndex(banknote$Status, banknote$kmeans)
mapping_kmeans <- bestMapping(banknote$Status, banknote$kmeans)
cm_kmeans <- mapping_kmeans$cm
acc_kmeans <- cm_kmeans$overall["Accuracy"]
sens_kmeans <- cm_kmeans$byClass["Sensitivity"]
spec_kmeans <- cm_kmeans$byClass["Specificity"]
distances_km <- t(apply(data_scaled, 1, function(x) {
apply(km_model$centers, 1, function(cent) sqrt(sum((x - cent)^2)))
}))
score_kmeans <- distances_km[,2] / rowSums(distances_km)
roc_obj_kmeans <- roc(banknote$Status, score_kmeans)
auc_kmeans <- auc(roc_obj_kmeans)
# K-medoids
ari_kmedoids <- adjustedRandIndex(banknote$Status, banknote$kmedoids)
mapping_kmedoids <- bestMapping(banknote$Status, banknote$kmedoids)
cm_kmedoids <- mapping_kmedoids$cm
acc_kmedoids <- cm_kmedoids$overall["Accuracy"]
sens_kmedoids <- cm_kmedoids$byClass["Sensitivity"]
spec_kmedoids <- cm_kmedoids$byClass["Specificity"]
medoids <- as.matrix(pam_model$medoids)
distances_kmed <- t(apply(data_scaled, 1, function(x) {
apply(medoids, 1, function(med) sqrt(sum((x - med)^2)))
}))
score_kmedoids <- distances_kmed[,2] / rowSums(distances_kmed)
roc_obj_kmedoids <- roc(banknote$Status, score_kmedoids)
auc_kmedoids <- auc(roc_obj_kmedoids)
# Gerarchico Agglomerativo
ari_hierAgg <- adjustedRandIndex(banknote$Status, banknote$hierarchical_agglomerative)
mapping_hierAgg <- bestMapping(banknote$Status, banknote$hierarchical_agglomerative)
cm_hierAgg <- mapping_hierAgg$cm
acc_hierAgg <- cm_hierAgg$overall["Accuracy"]
sens_hierAgg <- cm_hierAgg$byClass["Sensitivity"]
spec_hierAgg <- cm_hierAgg$byClass["Specificity"]
auc_hierAgg <- NA
# Gerarchico Divisivo
ari_hierDiv <- adjustedRandIndex(banknote$Status, banknote$hierarchical_divisive)
mapping_hierDiv <- bestMapping(banknote$Status, banknote$hierarchical_divisive)
cm_hierDiv <- mapping_hierDiv$cm
acc_hierDiv <- cm_hierDiv$overall["Accuracy"]
sens_hierDiv <- cm_hierDiv$byClass["Sensitivity"]
spec_hierDiv <- cm_hierDiv$byClass["Specificity"]
auc_hierDiv <- NA
# GMM
ari_gmm <- adjustedRandIndex(banknote$Status, banknote$gmm)
mapping_gmm <- bestMapping(banknote$Status, banknote$gmm)
cm_gmm <- mapping_gmm$cm
acc_gmm <- cm_gmm$overall["Accuracy"]
sens_gmm <- cm_gmm$byClass["Sensitivity"]
spec_gmm <- cm_gmm$byClass["Specificity"]
score_gmm <- gmm_model$z[,2]
roc_obj_gmm <- roc(banknote$Status, score_gmm)
auc_gmm <- auc(roc_obj_gmm)Infine, raccogliamo i risultati in un dataframe riassuntivo.
results_df <- data.frame(
Metodo = c("K-means", "K-medoids", "Gerarchico Agglomerativo", "Gerarchico Divisivo", "GMM"),
ARI = c(ari_kmeans, ari_kmedoids, ari_hierAgg, ari_hierDiv, ari_gmm),
Accuracy = c(acc_kmeans, acc_kmedoids, acc_hierAgg, acc_hierDiv, acc_gmm),
Sensitivity = c(sens_kmeans, sens_kmedoids, sens_hierAgg, sens_hierDiv, sens_gmm),
Specificity = c(spec_kmeans, spec_kmedoids, spec_hierAgg, spec_hierDiv, spec_gmm),
AUC = c(auc_kmeans, auc_kmedoids, auc_hierAgg, auc_hierDiv, auc_gmm)
)
print(results_df)## Metodo ARI Accuracy Sensitivity Specificity AUC
## 1 K-means 0.8456292 0.960 1.00 0.92 0.9973
## 2 K-medoids 0.9406018 0.985 1.00 0.97 0.9991
## 3 Gerarchico Agglomerativo 0.9602001 0.990 1.00 0.98 NA
## 4 Gerarchico Divisivo 0.9406015 0.985 0.99 0.98 NA
## 5 GMM 0.9799995 0.995 1.00 0.99 0.9999I risultati mostrano che il metodo GMM ha la migliore performance complessiva, con il più alto valore di ARI (0.98) e un’accuratezza del 99.5%. Il clustering gerarchico agglomerativo ha anch’esso buone prestazioni, mentre K-means e K-medoids mostrano una buona affidabilità ma con AUC leggermente inferiori. Analizzando la sensibilità e la specificità, si nota che la sensibilità – ovvero la capacità di identificare correttamente le banconote contraffatte – è generalmente molto alta per tutti gli algoritmi, indicando un’elevata capacità di individuare i falsi. Tuttavia, la specificità, ossia la capacità di riconoscere correttamente le banconote genuine, varia maggiormente tra i metodi. In particolare, K-means mostra una specificità inferiore (0.92) rispetto a K-medoids (0.97) e ai metodi gerarchici (0.98), suggerendo una maggiore tendenza a classificare erroneamente alcune banconote genuine come contraffatte; questo errore è meno grave di classificare una banconota contraffatta come vera.
3.1.5.3.2. Dataset ridotto
Nel seguente codice, applichiamo i diversi algoritmi di clustering al dataset ridotto (con le sole variabili ‘Diagonal’ e ‘Bottom’) con i migliori parametri calcolati in precedenza (fissando il numero di cluster a 2) e aggiungiamo al dataframe le colonne corrispondenti alle assegnazioni di ciascun metodo.
library(cluster)
data("banknote")
data <- banknote[, c("Bottom", "Diagonal")]
data_scaled <- scale(data)
set.seed(123)
km_model <- kmeans(data_scaled, centers = 2, nstart = 25)
pam_model <- pam(data_scaled, k = 2, nstart = 25)
#TODO: verificare se ward.D2 è il migliore
hc_model <- hclust(dist(data_scaled, method = "euclidean"), method = "ward.D2")
diana_model <- diana(data_scaled)
gmm_model <- Mclust(data_scaled, G = 2)
banknote$kmeans <- km_model$cluster
banknote$kmedoids <- pam_model$clustering
banknote$hierarchical_agglomerative <- cutree(hc_model, k = 2)
banknote$hierarchical_divisive <- cutree(as.hclust(diana_model), k = 2)
banknote$gmm <- gmm_model$classification
head(banknote)## Status Length Left Right Bottom Top Diagonal kmeans kmedoids
## 1 genuine 214.8 131.0 131.1 9.0 9.7 141.0 2 1
## 2 genuine 214.6 129.7 129.7 8.1 9.5 141.7 2 1
## 3 genuine 214.8 129.7 129.7 8.7 9.6 142.2 2 1
## 4 genuine 214.8 129.7 129.6 7.5 10.4 142.0 2 1
## 5 genuine 215.0 129.6 129.7 10.4 7.7 141.8 2 1
## 6 genuine 215.7 130.8 130.5 9.0 10.1 141.4 2 1
## hierarchical_agglomerative hierarchical_divisive gmm
## 1 1 1 1
## 2 1 1 1
## 3 1 1 1
## 4 1 1 1
## 5 1 1 1
## 6 1 1 1Nel blocco successivo, calcoliamo le metriche di valutazione come nel caso precedente.
# K-means
ari_kmeans <- adjustedRandIndex(banknote$Status, banknote$kmeans)
mapping_kmeans <- bestMapping(banknote$Status, banknote$kmeans)
cm_kmeans <- mapping_kmeans$cm
acc_kmeans <- cm_kmeans$overall["Accuracy"]
sens_kmeans <- cm_kmeans$byClass["Sensitivity"]
spec_kmeans <- cm_kmeans$byClass["Specificity"]
distances_km <- t(apply(data_scaled, 1, function(x) {
apply(km_model$centers, 1, function(cent) sqrt(sum((x - cent)^2)))
}))
score_kmeans <- distances_km[,2] / rowSums(distances_km)
roc_obj_kmeans <- roc(banknote$Status, score_kmeans)
auc_kmeans <- auc(roc_obj_kmeans)
# K-medoids
ari_kmedoids <- adjustedRandIndex(banknote$Status, banknote$kmedoids)
mapping_kmedoids <- bestMapping(banknote$Status, banknote$kmedoids)
cm_kmedoids <- mapping_kmedoids$cm
acc_kmedoids <- cm_kmedoids$overall["Accuracy"]
sens_kmedoids <- cm_kmedoids$byClass["Sensitivity"]
spec_kmedoids <- cm_kmedoids$byClass["Specificity"]
medoids <- as.matrix(pam_model$medoids)
distances_kmed <- t(apply(data_scaled, 1, function(x) {
apply(medoids, 1, function(med) sqrt(sum((x - med)^2)))
}))
score_kmedoids <- distances_kmed[,2] / rowSums(distances_kmed)
roc_obj_kmedoids <- roc(banknote$Status, score_kmedoids)
auc_kmedoids <- auc(roc_obj_kmedoids)
# Gerarchico Agglomerativo
ari_hierAgg <- adjustedRandIndex(banknote$Status, banknote$hierarchical_agglomerative)
mapping_hierAgg <- bestMapping(banknote$Status, banknote$hierarchical_agglomerative)
cm_hierAgg <- mapping_hierAgg$cm
acc_hierAgg <- cm_hierAgg$overall["Accuracy"]
sens_hierAgg <- cm_hierAgg$byClass["Sensitivity"]
spec_hierAgg <- cm_hierAgg$byClass["Specificity"]
auc_hierAgg <- NA
# Gerarchico Divisivo
ari_hierDiv <- adjustedRandIndex(banknote$Status, banknote$hierarchical_divisive)
mapping_hierDiv <- bestMapping(banknote$Status, banknote$hierarchical_divisive)
cm_hierDiv <- mapping_hierDiv$cm
acc_hierDiv <- cm_hierDiv$overall["Accuracy"]
sens_hierDiv <- cm_hierDiv$byClass["Sensitivity"]
spec_hierDiv <- cm_hierDiv$byClass["Specificity"]
auc_hierDiv <- NA
# GMM
ari_gmm <- adjustedRandIndex(banknote$Status, banknote$gmm)
mapping_gmm <- bestMapping(banknote$Status, banknote$gmm)
cm_gmm <- mapping_gmm$cm
acc_gmm <- cm_gmm$overall["Accuracy"]
sens_gmm <- cm_gmm$byClass["Sensitivity"]
spec_gmm <- cm_gmm$byClass["Specificity"]
score_gmm <- gmm_model$z[,2]
roc_obj_gmm <- roc(banknote$Status, score_gmm)
auc_gmm <- auc(roc_obj_gmm)Infine, raccogliamo i risultati in un dataframe riassuntivo.
results_df <- data.frame(
Metodo = c("K-means", "K-medoids", "Gerarchico Agglomerativo", "Gerarchico Divisivo", "GMM"),
ARI = c(ari_kmeans, ari_kmedoids, ari_hierAgg, ari_hierDiv, ari_gmm),
Accuracy = c(acc_kmeans, acc_kmedoids, acc_hierAgg, acc_hierDiv, acc_gmm),
Sensitivity = c(sens_kmeans, sens_kmedoids, sens_hierAgg, sens_hierDiv, sens_gmm),
Specificity = c(spec_kmeans, spec_kmedoids, spec_hierAgg, spec_hierDiv, spec_gmm),
AUC = c(auc_kmeans, auc_kmedoids, auc_hierAgg, auc_hierDiv, auc_gmm)
)
print(results_df)## Metodo ARI Accuracy Sensitivity Specificity AUC
## 1 K-means 0.9799995 0.995 0.99 1.00 0.9999
## 2 K-medoids 0.9799995 0.995 0.99 1.00 0.9999
## 3 Gerarchico Agglomerativo 0.9602000 0.990 0.99 0.99 NA
## 4 Gerarchico Divisivo 0.9406015 0.985 0.98 0.99 NA
## 5 GMM 0.9799995 0.995 1.00 0.99 0.9995Questa configurazione semplifica il modello senza perdere
accuratezza, dimostrando che ‘Bottom’ e ‘Diagonal’ sono le variabili
chiave per distinguere le classi:
- K-means, K-medoids e GMM ottengono ARI 0.98, accuracy 99.5% e AUC ~1,
confermando che queste due variabili sono altamente
discriminative.
- GMM raggiunge sensibilità 100%, identificando tutte le banconote
contraffatte, mentre K-means e K-medoids garantiscono specificità 100%,
evitando falsi positivi.
3.1.5.3.3. Dataset trasformato
Nel seguente codice, applichiamo i diversi algoritmi di clustering al dataset trasformato (con le variabili trasformate mediante tecnica PCA) con i migliori parametri calcolati in precedenza (fissando il numero di cluster a 2) e aggiungiamo al dataframe le colonne corrispondenti alle assegnazioni di ciascun metodo.
data <- banknote_pca[, -1]
data_scaled <- scale(data)
set.seed(123)
km_model <- kmeans(data_scaled, centers = 2, nstart = 25)
pam_model <- pam(data_scaled, k = 2, nstart = 25)
#TODO: verificare se ward.D2 è il migliore
hc_model <- hclust(dist(data_scaled, method = "euclidean"), method = "ward.D2")
diana_model <- diana(data_scaled)
gmm_model <- Mclust(data_scaled, G = 2)
banknote$kmeans <- km_model$cluster
banknote$kmedoids <- pam_model$clustering
banknote$hierarchical_agglomerative <- cutree(hc_model, k = 2)
banknote$hierarchical_divisive <- cutree(as.hclust(diana_model), k = 2)
banknote$gmm <- gmm_model$classification
head(banknote)## Status Length Left Right Bottom Top Diagonal kmeans kmedoids
## 1 genuine 214.8 131.0 131.1 9.0 9.7 141.0 2 1
## 2 genuine 214.6 129.7 129.7 8.1 9.5 141.7 1 2
## 3 genuine 214.8 129.7 129.7 8.7 9.6 142.2 1 2
## 4 genuine 214.8 129.7 129.6 7.5 10.4 142.0 1 2
## 5 genuine 215.0 129.6 129.7 10.4 7.7 141.8 1 2
## 6 genuine 215.7 130.8 130.5 9.0 10.1 141.4 1 1
## hierarchical_agglomerative hierarchical_divisive gmm
## 1 1 1 1
## 2 2 2 2
## 3 2 2 2
## 4 2 2 2
## 5 2 2 2
## 6 1 2 2Nel blocco successivo, calcoliamo le metriche di valutazione come nei casi precedenti.
# K-means
ari_kmeans <- adjustedRandIndex(banknote$Status, banknote$kmeans)
mapping_kmeans <- bestMapping(banknote$Status, banknote$kmeans)
cm_kmeans <- mapping_kmeans$cm
acc_kmeans <- cm_kmeans$overall["Accuracy"]
sens_kmeans <- cm_kmeans$byClass["Sensitivity"]
spec_kmeans <- cm_kmeans$byClass["Specificity"]
distances_km <- t(apply(data_scaled, 1, function(x) {
apply(km_model$centers, 1, function(cent) sqrt(sum((x - cent)^2)))
}))
score_kmeans <- distances_km[,2] / rowSums(distances_km)
roc_obj_kmeans <- roc(banknote$Status, score_kmeans)
auc_kmeans <- auc(roc_obj_kmeans)
# K-medoids
ari_kmedoids <- adjustedRandIndex(banknote$Status, banknote$kmedoids)
mapping_kmedoids <- bestMapping(banknote$Status, banknote$kmedoids)
cm_kmedoids <- mapping_kmedoids$cm
acc_kmedoids <- cm_kmedoids$overall["Accuracy"]
sens_kmedoids <- cm_kmedoids$byClass["Sensitivity"]
spec_kmedoids <- cm_kmedoids$byClass["Specificity"]
medoids <- as.matrix(pam_model$medoids)
distances_kmed <- t(apply(data_scaled, 1, function(x) {
apply(medoids, 1, function(med) sqrt(sum((x - med)^2)))
}))
score_kmedoids <- distances_kmed[,2] / rowSums(distances_kmed)
roc_obj_kmedoids <- roc(banknote$Status, score_kmedoids)
auc_kmedoids <- auc(roc_obj_kmedoids)
# Gerarchico Agglomerativo
ari_hierAgg <- adjustedRandIndex(banknote$Status, banknote$hierarchical_agglomerative)
mapping_hierAgg <- bestMapping(banknote$Status, banknote$hierarchical_agglomerative)
cm_hierAgg <- mapping_hierAgg$cm
acc_hierAgg <- cm_hierAgg$overall["Accuracy"]
sens_hierAgg <- cm_hierAgg$byClass["Sensitivity"]
spec_hierAgg <- cm_hierAgg$byClass["Specificity"]
auc_hierAgg <- NA
# Gerarchico Divisivo
ari_hierDiv <- adjustedRandIndex(banknote$Status, banknote$hierarchical_divisive)
mapping_hierDiv <- bestMapping(banknote$Status, banknote$hierarchical_divisive)
cm_hierDiv <- mapping_hierDiv$cm
acc_hierDiv <- cm_hierDiv$overall["Accuracy"]
sens_hierDiv <- cm_hierDiv$byClass["Sensitivity"]
spec_hierDiv <- cm_hierDiv$byClass["Specificity"]
auc_hierDiv <- NA
# GMM
ari_gmm <- adjustedRandIndex(banknote$Status, banknote$gmm)
mapping_gmm <- bestMapping(banknote$Status, banknote$gmm)
cm_gmm <- mapping_gmm$cm
acc_gmm <- cm_gmm$overall["Accuracy"]
sens_gmm <- cm_gmm$byClass["Sensitivity"]
spec_gmm <- cm_gmm$byClass["Specificity"]
score_gmm <- gmm_model$z[,2]
roc_obj_gmm <- roc(banknote$Status, score_gmm)
auc_gmm <- auc(roc_obj_gmm)Infine, raccogliamo i risultati in un dataframe riassuntivo.
results_df <- data.frame(
Metodo = c("K-means", "K-medoids", "Gerarchico Agglomerativo", "Gerarchico Divisivo", "GMM"),
ARI = c(ari_kmeans, ari_kmedoids, ari_hierAgg, ari_hierDiv, ari_gmm),
Accuracy = c(acc_kmeans, acc_kmedoids, acc_hierAgg, acc_hierDiv, acc_gmm),
Sensitivity = c(sens_kmeans, sens_kmedoids, sens_hierAgg, sens_hierDiv, sens_gmm),
Specificity = c(spec_kmeans, spec_kmedoids, spec_hierAgg, spec_hierDiv, spec_gmm),
AUC = c(auc_kmeans, auc_kmedoids, auc_hierAgg, auc_hierDiv, auc_gmm)
)
print(results_df)## Metodo ARI Accuracy Sensitivity Specificity AUC
## 1 K-means 0.8830121 0.970 0.97 0.97 0.9987
## 2 K-medoids 0.8272389 0.955 1.00 0.91 0.9946
## 3 Gerarchico Agglomerativo 0.5160469 0.860 0.95 0.77 NA
## 4 Gerarchico Divisivo 0.9212040 0.980 0.98 0.98 NA
## 5 GMM 0.9212042 0.980 0.99 0.97 0.9983Rispetto ai risultati con il dataset completo o ridotto, l’uso di solo 3 variabili ottenute con PCA ha ridotto la performance globale, in particolare per K-medoids e Gerarchico Agglomerativo. Sebbene l’accuratezza resti elevata, l’ARI e la specificità mostrano una leggera diminuzione, il che suggerisce che la riduzione delle variabili ha compromesso la capacità di questi algoritmi di catturare la struttura dei dati. In particolare, l’algoritmo Gerarchico Agglomerativo mostra una performance molto più bassa in termini di ARI e accuratezza, indicando che l’algoritmo è più sensibile alla perdita di variabili discriminanti.
In conclusione, dal punto di vista della misurazione delle performance con tecniche esterne (avendo a disposizione la vera classe delle osservazioni), l’algoritmo GMM è risultato il migliore nei tre diversi dataset, anche se, tutto sommato, gli altri algoritmi non hanno ottenuto dei cattivi risultati.